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Abstract 

We performed a comprehensive study of the spike autosolitons: self-sustained soli- 
tary inhomogeneous states, in the classical reaction-diffusion system — the Gray- 
Scott model. We developed singular perturbation techniques based on the strong 
separation of the length scales to construct asymptotically the solutions in the form 
of a one-dimensional static autosolitons, higher-dimensional radially-symmetric static 
autosolitons, and two types of traveling autosolitons. We studied the stability of the 
static autosolitons in one and three dimensions and analyzed the properties of the 
static and the traveling autosolitons. 
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1 Introduction 



Self-organization and pattern formation in nonequilibrium systems are among 
the most fascinating phenomena in nonlinear physics [1-11]. Pattern formation 
is observed in various physical systems including aero- and hydrodynamic sys- 
tems; gas and electron-hole plasmas; various semiconductor, superconductor 
and gas-discharge structures; some ferroelectric, magnetic and optical media; 
combustion systems (see, for example, [5,9-14]), as well as in many chemical 
and biological systems (see, for example, [1-7,15]). 

Self-organization is often associated with the destabilization of the uniform 
state of the system [1,2,5,10,11]. At the same time, when the uniform state 
of the system is stable, by applying a sufficiently strong perturbation one can 
excite large-amplitude patterns, including autosolitons (ASs) — self-sustained 
solitary inhomogeneous states [8-11,16-19]. Autosolitons are the elementary 
objects in open dissipative systems away from equilibrium. They share the 
properties of both solitons and traveling waves (or autowaves, as they are also 
referred to [2,6]). They are similar to solitons since they are localized objects 
whose existence is due to the nonlinearities of the system. On the other hand, 
from the physical point of view they are essentially different from solitons 
in that they are dissipative structures, that is, they are self-sustained objects 
which form in strongly dissipative systems as a result of the balance between 
the dissipation and pumping of energy or matter. This is the reason why, 
in contrast to solitons, their properties are independent of the initial condi- 
tions and are determined primarily by the nonlinearities of the system [8-11]. 
ASs can be static, pulsating, or traveling. As a result of their various insta- 
bilities, these simplest localized patterns can spontaneously transform into 
complex space-filling static or dynamic patterns, including complex pulsat- 
ing and traveling patterns, or spatio-temporal chaos [8-14,18-33]. Thus, it is 
the destabilization of the ASs that is the main source of self-organization in 
nonequilibrium systems with the stable homogeneous state. 

Real physical, chemical, and biological systems exhibiting pattern formation 
and self-organization are extremely complicated, so simplified models are used 
to describe these phenomena. A prototype model of this kind is a pair of 
reaction-diffusion equations of the activator-inhibitor type 



where 9 is the activator, r\ is the inhibitor, r e , I and t v , L are the time and the 
length scales of the activator and the inhibitor, respectively; A is the control 
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l 2 A9-q{6^A), 



(1.1) 




L 2 Ar)-Q(9,r],A), 
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(bifurcation) parameter; q and Q are certain nonlinear functions representing 
the activation and the inhibition processes. Examples of these equations for 
various physical systems are given in [9-13,25] where the physical meaning 
of the variables 9 and r\ and the nature of the activation and the inhibition 
processes are discussed. The well-known Brusselator [1] and the Gray-Scott 
[34] models of autocatalytic chemical reactions, the classical Gierer-Meinhardt 
model of morphogenesis [35], the FitzHugh-Nagumo [36] and the piecewise- 
linear Rinzel-Keller model [37] for the propagation of pulses in the nerve fibers 
are all special cases of Eqs. (1.1) and (1.2). 

The fact that 9 is the activator means that for certain parameters the uniform 
fluctuations of 9 will grow when the value of 77 is fixed. From the mathematical 
point of view, this is given by the condition [9-11,18] 



for certain values of 9 and 77. On the other hand, the fact that r? is the inhibitor 
means that its own fluctuations decay and that it damps the fluctuations of 
the activator. Mathematically, these conditions are expressed by [9-11,18] 



for all values of 9 and rj, provided that the derivatives in Eq. (1.4) do not 
change sign. 

Kerner and Osipov showed [8-11,16-18,25] that the properties of the patterns 
and the self-organization scenarios in systems described by Eqs. (1.1) and 
(1.2) are chiefly determined by the parameters e = l/L and a = Tg/r v and the 
shape of the nullcline of the equation for the activator, that is, the dependence 
f](9) given by the equation q(9,rj,A) = for A = const. They demonstrated 
that depending on the shape of the activator nullcline all systems involved 
can be divided into two fundamentally different classes: N-systems, for which 
the nullcline is N- or inverted N-shaped and, A- or V-systems, for which the 
nullcline is A- or V-shaped, respectively (see Fig. 1). 

Most works devoted to the description of pattern formation on the basis of Eqs. 
(1.1) and (1.2) deal with N-systems. In N-systems the equation q(9,rj,A) = 
has three roots: 9i,9 2 , and 9 3 , for given values of A and 77. The roots 9\ and 
#3 correspond to the stable states and 9 2 corresponds to the unstable state 
in the system with rj = const. It is easy to see that the FitzHugh-Nagumo 
and the piecewise-linear models belong to N-systems. For these models it was 
proved [27,37] that Eqs. (1.1) and (1.2) with L = and a = tq/t^ <C 1 have 
solutions in the form of the traveling waves (also called autowaves [2,6], or 
traveling ASs [8-11]). In [16-19] it was shown that in another limit L 3> I 



q' e <0 



(1.3) 



Q;>0, q' v Q' e <0 



(1.4) 
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(or, more precisely, when e = I/L < 1 and a > 1) Eqs. (1.1) and (1.2) admit 
solutions in the form of the stable static patterns including ASs (see also [9- 
11]). Furthermore, it was shown that in systems with e < 1 and a C 1 one 
can excite static, pulsating, and traveling patterns [8-11,18,25,28-33]. The 
characteristic velocity of the traveling patterns in N-systems does not exceed 
the value of order l/r e [9,11,18,27]. 

It is important to emphasize that e or a are the natural small parameters in 
these system. Their smallness is in fact a necessary condition for the feasibility 
of any patterns [9-11]. Indeed, if the inverse were true, that is, if both the 
characteristic time and length scales of the variation of the inhibitor were 
much smaller than those of the activator, the inhibitor would easily damp all 
the deviations of the activator from the homogeneous steady state, making the 
formation of any kinds of persistent patterns impossible. On the other hand, 
the fact that we must have either e < 1 or a < 1 implies that it is advantageous 
to consider the extreme cases of e 1 or a <C 1. These conditions, in turn, will 
result in a significant simplification of the original highly nonlinear problem. 
Recently, this kind of approach has been successfully applied to a variety of 
problems (see, for example, [38]). 

In N-systems with e < 1 the static ASs and other patterns are essentially the 
domains of high and low values of the activator separated by the interfaces 
(walls) where 9 varies sharply over a distance of order I from one stable state 
Q\ to the other 63. The characteristic size C s of these domain patterns lies in 
the range / <C C s < L and their amplitude (the value 9\— 63) is determined 
by the form of q and Q and becomes independent of e as e — > [9-11,17-22]. 
The properties of these patterns are essentially determined by the dynamics 
of their interfaces and the interaction between them. So, the majority of the 
theoretical work devoted to the description of the complex domain patterns 
in N-systems in fact developed an interfacial dynamics approach [9-11,17- 
25,29-33,39-41]. On the basis of this approach, Kerner and Osipov developed 
a theory of instabilities of the domain patterns in the general N-systems in 
one dimension [8-11]. More recently, we extended this theory for arbitrary do- 
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main patterns in higher-dimensional systems and extensively studied various 
pattern formation scenarios in these systems [21,22,40,41]. 

At the same time, there are many physical, chemical and biological systems 
for which the activator nullcline is A- or V-shaped [Fig. 1(b)]. In this case the 
equation q(6, i],A) = for given A and rj has only two roots: Q\ corresponding 
to the stable state, and 82 corresponding to the unstable state in the sys- 
tem with rj = const [9-11,16]. Among A-systems are many semiconductor and 
gas discharge structures, electron-hole and gas plasmas, radiation heated gas 
mixtures (see, for example, [9-12,16,25]). It is not difficult to see that the Brus- 
selator and the Gray-Scott models are A-systems, and the Gierer-Meinhardt 
model is a V-system. 

Kerner and Osipov qualitatively showed that in A-systems the so-called spike 
ASs and more complex spike patterns can be excited [9,11,16,42,43]. They 
were the first to analyze the static spike ASs and strata in the Brusselator, 
the Gierer-Meinhardt model, and the electron-hole plasma [16,42]. They found 
that when e C 1 and a > 1, the one- dimensional static spike AS can have 
small size of order / and huge amplitude which goes to infinity as e — > 0. 
Dubitskii, Kerner, and Osipov formulated the asymptotic procedure for finding 
the stationary solutions in A-systems for sufficiently small e [11,42]. Recently, 
we showed that in another limiting case a C 1 and e ^> 1 one can excite 
the one-dimensional traveling spike AS which also has small size and whose 
amplitude goes to infinity as a — > [44]. We also showed that, in contrast 
to the traveling patterns in N-systems, the velocity of this one-dimensional 
traveling spike AS can have huge values (c ^> 1/tq) and that the inhibitor 
distribution varies stepwise in the front of the spike. Thus, one can see that 
the properties of the spike patterns forming in A-systems differ fundamentally 
from those of the domain patterns forming in N-systems. In particular, since 
the interface connecting the two stable states at rj = const does not exist 
in A- systems, the size of the spike should be of the order of the smallest 
system length scale. For this reason, the concept of the interfacial dynamics 
developed for the domain patterns in N-systems is generally inapplicable to the 
description of spike patterns. Some properties of the one-dimensional static 
spike ASs and the main types of their instabilities in the simplified version 
of the Gray-Scott model have recently been studied by Osipov and Severtsev 
[45]. 

Spike patterns including the spike ASs are observed experimentally in the 
nerve tissue [46], chemical reactions [5,47], electron-hole plasma [48], gas- 
discharge structures [49] , as well as numerically in the simulations of the Brus- 
selator, the Gierer-Meinhardt, and the Gray-Scott models [1,15,26,35,50]. At 
the same time, there is a only a limited number of theoretical studies of these 
patterns. Moreover, many aspects of these patterns including the stability of 
the spike ASs and their properties in higher dimensions as well as the spon- 
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taneous transitions between the static, the pulsating and the traveling spike 
ASs have not been studied at all. So far, there have been no general methods 
for dealing with spike patterns. 

In the present paper we develop asymptotic methods for the description of 
the spike patterns and study their major properties in arbitrary dimensions 
as well as their shape and stability. We find the conditions of the spontaneous 
transitions between different types of the spike ASs and study the scenarios of 
the formation of the spike ASs and more complex spike patterns in one- and 
two-dimensional systems. To be specific, we consider the Gray-Scott model of 
an autocatalytic chemical reaction, which possesses a number of advantages. 
First, it is one of the rarest models for which in many cases one can obtain ex- 
act results. Second, a lot of the numerical studies of this model were performed 
recently [15,26,50,51]. Finally, the Gray-Scott model possesses a particularly 
simple set of nonlinearities, so one can expect a certain degree of universality 
in the pattern formation scenarios exhibited by it. 

The outline of our paper is as follows. In Sec. 2 we introduce the model we will 
study, in Sec. 3 we asymptotically construct the one- dimensional static AS, 
and the two- and the three-dimensional radially-symmetric ASs, in Sec. 4.1 
we asymptotically construct the solutions in the form of the two types of trav- 
eling spike ASs, in Sec. 5 we analyze the stability of the one-dimensional and 
the higher-dimensional radially-symmetric static ASs and show the existence 
of various instabilities, in Sec. 6 we compare our results with the numerical 
simulations of the one-dimensional system, and in Sec. 7 we do that for the 
two-dimensional system, in Sec. 8 we discuss the works of other authors on the 
Gray-Scott model in light of our results, and in Sec. 9 we give the summary 
of our work and draw conclusions. 



2 The model 

The Gray-Scott model describes the kinetics of a simple autocatalytic reaction 
in an unstirred flow reactor. The reactor is a narrow space between two porous 
walls. Substance Y whose concentration is kept fixed outside of the reactor is 
supplied through the walls into the reactor with the rate ko and the products 
of the reaction are removed from the reactor with the same rate. Inside the 
reactor Y undergoes the reaction involving an intermediate species X: 



2X + Y -^ 3X, 
X % inert. 



(2.1) 
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The first reaction is a cubic autocatalytic reaction resulting in the self-production 
of species X; therefore, X is the activator species. On the other hand, the 
production of X is controlled by species Y, so Y is the inhibitor species. The 
equations of chemical kinetics which describe the spatiotemporal variations of 
the concentrations of X and Y in the reactor and take into account the supply 
and the removal of the substances through the porous walls take the following 
form [34]: 

8X 

— = -(k + k 2 )X + k 1 X 2 Y + D x AX, (2.3) 
at 

BY 

— = k (Y — Y) — hX 2 Y + D Y AY, (2.4) 
at 

where now X and Y are the concentrations of the activator and the inhibitor 
species, respectively, Y is the concentration of Y in the reservoir, A is the 
two-dimensional Laplacian, and D x and D Y are the diffusion coefficients of 
X and Y. 

In order to be able to understand various pattern formation phenomena in a 
system of this kind, it is crucial to introduce the variables and the time and 
length scales that truly represent the physical processes acting in the system. 
The first and the most important is the choice of the characteristic time scales. 
These are primarily dictated by the time constants of the dissipation processes. 
For Y this is the supply and the removal with the rate ho, whereas for X this 
is the removal from the system and the decay via the second reaction with the 
total rate ko + k 2 . The natural way to introduce the dimensionless inhibitor 
concentration is to scale it with Y . Since we want to fix the time scale of the 
variation of the inhibitor (with the fixed activator), we will rescale X in such 
a way that the reaction term in Eq. (2.4) will generate the same time scale as 
the dissipative term. This leads to the following dimensionless quantities: 



(k \ 1/2 

e = X/X , r ] = Y/Y , ^o=(^J • (2.5) 
The characteristic time and length scales for these quantities are 



re = (k + k 2 ) \ r v = k 1 , (2.6) 
l = (D x r e ) 1/2 , L = (Dyr,) 1 ' 2 . (2.7) 

Naturally, one should require the positivity of 6 and rj. 

If we now write Eqs. (2.3) and (2.4) in the dimensionless form, we will arrive 
at the following set of equations: 
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Fig. 2. The nullclines of Eqs. (2.8) and (2.9) for A = 1 (a) and A = 3 (b). 
f)9 

r e — = l 2 A9 + A9 2 r ] -9, (2.8) 
at 

T v ^ = L 2 Ar)-9 2 r) + l-r), (2.9) 
where we introduced a dimensionless parameter 



Y h l/2 h 112 

A= tp ?. . (2.10) 

One can see from Eqs. (2.8) and (2.9) that tq and r v are in fact the character- 
istic time scales, and / and L the characteristic length scales of the variation 
of small deviations of 9 and i] from the stationary homogeneous state 9 = 9h 
and 7] = r] h : 



9 h = 0, Vh = l. (2.11) 

Thus, the system is characterized by only three dimensionless parameters: 
a = Tfl/rjj, e = l/L, and A. As can be seen from Eq. (2.8), the parameter A is 
the dimensionless strength of the activation process, that is, it describes the 
degree of deviation of the system from thermal equilibrium. With all this, Eqs. 
(2.8) and (2.9) are reduced to the form of Eqs. (1.1) and (1.2). Notice that 
the system given by Eqs. (2.8) and (2.9) is indeed a system of the activator- 
inhibitor type: the condition in Eq. (1.3) is satisfied for 9 > an d the 
conditions in Eq. (1.4) are satisfied with q' n < and Q' e > for all 9 > and 
7] > 0. 

The nullclines of Eqs. (2.8) and (2.9) are shown in Fig. 2. From this figure one 
can see that the nullcline of the equation for the activator has degenerate A- 
form. It consists of two separate branches: 9 = and 9 — One can easily 
check that for < A < 2 there is only one stationary homogeneous state 
given by Eq. (2.11), whereas for A > 2 two extra stationary homogeneous 
states exist 
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A T V^A A±y/A?=Z 
Vh2,3 ^ ' ^2,3 ~~ 2A ' \ lAl ) 

The stability analysis of these homogeneous states shows that for e <C 1 or 
a <C 1 the homogeneous state 9 = 9h2, V = Vh2 is always unstable. For 
f < 1 the homogeneous state 9 = 9^3, rj = 77^3 is unstable with respect to the 
Turing instability if A < 0.41e _1 . For a 1 it is unstable with respect to the 
homogeneous oscillations (Hopf bifurcation) if OAla' 1 ^ 2 < A < a' 1 / 2 , or it is 
an unstable node if A < 0.41a -1 / 2 . On the other hand, the homogeneous state 
9 = Oh, rj = 7]^ is stable for all values of the system's parameters. The latter is 
simple to understand: in order for the reaction to begin there has to be at least 
some amount of the activator put in at the start. Equivalently, the fact that 
the homogeneous state in Eq. (2.11) is stable for all values of the parameter 
A (for an arbitrary deviation from thermal equilibrium) is the consequence of 
the degeneracy of the nullcline of Eq. (2.8). Thus, self-organization associated 
with the Turing instability of the homogeneous state 9^ = and 77^ = 1 is not 
realized in the Gray-Scott model. In such a stable homogeneous system any 
inhomogeneous pattern, including the ASs, can only be excited by a sufficiently 
strong localized stimulus. In turn, self-organization will occur as a result of 
the instabilities of the large-amplitude patterns already present in the system. 

Note that in the opposite case e> 1 and a»l the dynamics of the system 
becomes dramatically simpler. Indeed, if we put both L and r v to zero, from 
Eq. (2.9) we get a local relationship rj = j^p- Substituting this back to Eq. 
(2.8), we obtain 



This equation possesses a simple variational structure 



T ,| = _^, r = J^(mr- M + A m »n t + *). (2,4) 

For A < 2 the functional T has a unique global minimum at 9 = 9h = 0, so 
any initial condition will relax to the homogeneous state Oh- For A > 2 there 
are two stable homogeneous states 9 = 9 h and 9 = 6^3 (see above), so it is 
possible to have the waves of switching from one homogeneous state to the 
other [5]. It is easily checked that for 2 < A < 2.18 the dominant homogeneous 
state is 9 h , while for A > 2.18 the dominant homogeneous state is 6^3. 

In the case of e 1 and a < 1 the largest length scale in the system is L and 
the longest time scale is r v , so it is natural to scale length and time with L 
and t v , respectively. In these units Eqs. (2.8) and (2.9) will take the following 
form: 
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m 

di 
drj 

di 



Ar] - d 2 rj + 1 - 77. 



e 2 A9 + A9 2 r] - 9, 



(2.15) 



(2.16) 



We will assume that the problem is defined on the sufficiently large domain 
with neutral boundary conditions. Notice that the kinetic model used to arrive 
at Eqs. (2.15) and (2.16) imposes a restriction a < 1 [see Eq. (2.6)]. Also, in 
the derivation we assumed that the system is essentially two-dimensional. For 
the sake of generality, in the following we will allow a to take arbitrary values 
and will work with the arbitrarily-dimensional Gray-Scott model. 



3 Static spike autosolitons 

Let us now study the simplest possible stationary pattern in the Gray-Scott 
model — the static spike AS. According to the general qualitative theory, 
these ASs form in A-systems when e C 1 [9-11,16]. The condition e < 1 will 
therefore be assumed throughout this section. 



3.1 One- dimensional static spike autosoliton 

We begin with the analysis of the one-dimensional static spike AS. In the 
Gray-Scott model it is described by the following equations 



Since e C 1, there is a strong separation of the length scales in the AS [9- 
11,16]. One can separate the spike region where the distribution of 9 varies 
on the length scale of e, and the periphery of the AS where 77 decays into the 
homogeneous state r\ h = 1 on the length scale of order 1. One can use this 
separation of the length scales to construct a singular perturbation theory 
which describes the distributions in the form of the static one-dimensional 
spike AS [42]. But before we do that, it is instructive to use a more qualitative 
approach which will give us an idea about the scaling of the main parameters 
of the AS and its qualitative shape. As will be seen below, this approach works 
when e 1 / 2 < A < 1. 




(3.1) 



(3.2) 
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3.1.1 Case A ~ e 1 / 2 : autosoliton collapse 



According to this approach [9-11,16], one assumes that the value of i] inside 
the spike (on the length scale of e) is close to a constant. This is a reasonable 
assumption as long as r\ ^> e in the spike since the characteristic length scale 
of the variation of r\ is 1. Let us denote this constant value of rj as i] s . Then, 
Eq. (3.1) with rj — r) s can be solved exactly. Its solution has the form 



f)(.r) sb 2 (£j with^ —. (3.3) 



On the other hand, the distribution of 6 given by Eq. (3.3) acts in Eq. (3.2) 
as a 5-function, so away from the spike the distribution of rj is given by 



iris) = 1 - e-W . (3-4) 

Now, matching this solution for r)(x) with the condition that 77(0) = r) s , we 
obtain the following expressions 



m — —77 



3A 



Vs = 



A 2 
2A 2 




(3.5) 



where 



A b = V12e. (3.6) 

Note that these results were also obtained in [51] by applying Melnikov analysis 
to Eqs (3.1) and (3.2). Similar results for the simplified version of the Gray- 
Scott model were obtained in [45]. 

From Eq. (3.6) one can see that at A < A b the solution in the form of the spike 
AS does not exist. When A > A b there are two solutions: the one corresponding 
to the plus sign has larger amplitude and the one corresponding to the minus 
sign has smaller amplitude. As was shown by Kerner and Osipov, the solutions 
that have smaller amplitude are always unstable [9-11], so the only interesting 
solution corresponds to the plus sign in Eq. (3.5). This solution is precisely 
the static spike AS. The numerical simulations of Eqs. (2.15) and (2.16) show 
that if the value of A is lowered, at A = A b a stable one-dimensional static 
spike AS collapses into the homogeneous state (see Sec. 6). 

Let us look more closely at the parameters of the static spike AS and the 
conditions of validity of the approximations made in the preceding paragraphs. 
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As can be seen from Eqs. (3.3) and (3.4), the distribution of the activator 
indeed has a form of the spike whose characteristic width is of order e, and 
the distribution of the inhibitor varies on the much larger length scale of order 
1. Also, according to Eq. (3.5), the amplitude of the spike at A close to A b is of 
order e -1 / 2 ^> 1 (and can in fact have huge values as e gets smaller) and grows 
as the value of A increases. At A ~ 1 the amplitude 9 m ~ e _1 . These features 
fundamentally differ the AS forming in A-systems from the AS in N-systems. 

Recall that in the derivation we neglected the variation of the inhibitor inside 
the spike. Since the characteristic length of the variation of r\ is of order 1, this 
means that the value of rj = r/ s in the center of the AS must be much greater 
than e. According to Eq. (3.5), this is indeed the long as i C 1, so 

the solution obtained above is a good approximation to the actual solution 
in this case. Also, in this case one can easily calculate the distribution of rj 
in the spike. To do this, we note that, according to Eq. (3.5), for A C 1 we 
have 9 2 r\ 3> 1 in the spike, so the last two terms in Eq. (3.2) can be neglected. 
Since the variation of rj in the spike is much smaller compared to rj s , we can 
put rj = rj s in the right-hand side of Eq. (3.2). Then, substituting 9 from Eq. 
(3.3) into this equation, after simple integration we obtain an expression for 
f] in the spike region 

, , 4e 2 9 2 n ri s / , x 1 , 9 x \ , 

r}{x) =r) s + 2 In cosh — + - tanh 2 — . 3.7 

3 V 2e 2 ZeJ 

3.1.2 Case A ~ 1: local breakdown 

On the other hand, according to Eqs. (3.5), when A ~ 1, we have 

#max ~ e~\ Tfcnin ~ 6, (3-8) 

and the approximation used by us ceases to be valid. However, it is clear that 
qualitatively the character of the solution should not change even for these 
values of A. Therefore, we can still assume that the spike of the AS has the 
width of order e and that the values of the activator and the inhibitor scale 
the same way as those in Eq. (3.8). With all this in mind, we are now able 
to introduce singular perturbation expansion and separate the "sharp" distri- 
butions (inner solutions) that vary on the length scale of e and the "smooth" 
distributions (outer solutions) that vary on the length scale of order 1. 

At distances much greater than e away from the spike (in the outer region) the 
value of 9 is exponentially zero,Q so the equation for the smooth distributions 

1 This follows from the fact that in the region of the smooth distributions 9 and rj 
are related locally through the equation q(9,rj) = and therefore must lie on the 
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becomes 



§ + 1-1 = 0, (3.9) 

with the boundary condition in the spike 77(0) = (to order e) and 77 = r]h at 
infinity. This immediately gives us the smooth distribution of 77 



rj( x ) = 1 - . (3.10) 



Let us scale the activator and the inhibitor according to Eq. (3.8) and intro- 
duce the stretched variable £: 



e = e 9, f) = e-\ i = -^. (3.11) 

Using these variables, after a little algebra we can write Eqs. (3.1) and (3.2) 

as 



% + A9 2 fi - 9 = 0, (3.12) 
ite = A- 1 {9-9x) 1 (3.13) 

where we kept only the leading terms. In this equation 9^ denotes the sec- 
ond derivative with respect to £. Thus, the solution of Eqs. (3.12) and (3.13) 
properly matched with the smooth distribution, given by Eq. (3.10), will give 
the sharp distributions of the activator and the inhibitor in the spike. 

The matching of the sharp and the smooth distributions is performed by noting 
that, according to Eq. (3.11), to order e we have rj = and rj x ~ 1 for e <C 
\x\ <C I. Therefore, it is the derivative of rj obtained from the sharp distribution 
at |£| ^> 1 that must coincide with that of the smooth distribution for \x\ 1. 
This condition is obtained by imposing the boundary condition ^(±00) = ±1 
in Eq. (3.13) [see Eq. (3.10)]. One can obtain an integral representation of 
this boundary condition by integrating Eq. (3.13) over £. Let us introduce the 
variables 



9 = j, v = f, + j. (3.14) 
Then, this integral condition takes the form 



stable branch of the nullcline of the equation for the activator [9—11], which in the 
case of the Gray-Scott model gives especially simple relation: 6 = 



14 



+oo 

y^ = 2|A|- 1/2 , (3.15) 

— oo 

where A is a constant that should be equal to —1 (the reason for introducing 
this coefficient will be explained in the following paragraph). In terms of the 
new variables Eq. (3.13) becomes especially simple 

% = |A|0, (3.16) 

where A is the same constant. Notice that Eq. (3.16) has an obvious symmetry 
which allows us to add an arbitrary constant to fj, so we can replace fj — > 
fj + fj s , where fj satisfies the condition 77(0) = %(0) = and thus is uniquely 
determined by 9, and fj s is an arbitrary constant. 

To analyze Eqs. (3.12) and (3.13), it is convenient to rewrite them as a non- 
linear eigenvalue problem 



d 2 



de , = * (3.17) 

V(S) = -A 2 e(f} a + f}-8), (3.18) 

where fj is in turn related to 9 via Eq. (3.16). Then, since 9 is positive for all 
£ and therefore has no nodes, the solution in the form of the static spike AS 
will correspond to the lowest bound state of the operator in Eq. (3.17) with 
A = — 1. The latter is achieved by adjusting the value of fj s . 

The nonlinear eigenvalue problem given by Eqs. (3.17) and (3.18) together 
with Eqs. (3.15) and (3.16) with fixed fj s can be solved iteratively. Indeed, 
for a given potential well V there is a unique eigenvalue A and a unique 
eigenfunction 9 (up to normalization) that correspond to the lowest bound 
state of the Schrodinger operator in Eq. (3.17). Equation (3.15) gives a unique 
normalization for 9, which then uniquely determines fj through Eq. (3.16). 
Knowing the distributions 9 and fj, one can then reconstruct the potential V , 
thus defining an iterative map. It is convenient to think of the solutions of the 
nonlinear eigenvalue problem as fixed points of this iterative map. 

Observe that the nonlinear eigenvalue problem is invariant with respect to the 
following transformation 



A^6 2 A, A — > bA, (3.19) 

where b is an arbitrary positive constant. It is clear that if one knows a solution 
of the nonlinear eigenvalue problem with certain A, one can obtain a solution 



15 




Fig. 3. Qualitative form of the potentials Vq and V\ from Eqs. (3.20). 

of Eqs. (3.12) and (3.13) by simply using the symmetry transformation in 
Eq. (3.19) with b = \\\' 1/2 , so there is in fact a one-to-one correspondence 
between the solutions of the nonlinear eigenvalue problem with arbitrary A 
and its solution with A = — 1 which corresponds to the sharp distributions. 

Since we are interested in the lowest bound state whose eigenvalue is equal to 
— 1, the characteristic length scale of the variation of 9 and, according to Eq. 
(3.16), of fj and V as well, is of order 1. Equation (3.15) fixes the normalization 
of 9, so we must have 9 ~ 1. Notice that in view of Eq. (3.16) and the fact 
that 9 ~ 1, we must always have fj ~ 1. 

Let us write the potential V in Eq. (3.17) and (3.18) as a sum of two parts: 
V — Vq + Vi, where 



V = -A 2 9fi s , V 1 = -A 2 9{fj - 9). (3.20) 

From the qualitative form of r/(£) [see, for example, Eq. (3.7)] it is easy to see 
that the potential Vq has the form of a simple potential well, while V 1 has the 
form of a double well (see Fig. 3). 

In order for the operator in Eq. (3.17) to have the lowest bound state with 
A = — 1 the potential V must have the depth of order 1. When A C 1, the 
function ??(£) must be chosen in such a way that it compensates for the small 
factor of A in Eq. (3.20). However, since fj ~ 1, this can only be achieved by 
choosing fj s ~ A~ 2 ^> 1. This means that we will have Vq ^> V\. If one neglects 
V\ compared to Vo, one can solve the nonlinear eigenvalue problem exactly. 
This solution will be 



^(0 = ^cosh- 2 (jj, ^ = -|. (3.21) 

The potential V\ can then be treated as a perturbation which will give cor- 
rections to 9 and fj, so one should not expect any qualitative changes in the 
behavior of the solution for A < 1. 
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Fig. 4. Dependence X(fj s ) for the one-dimensional static AS. Results of the numerical 
solution of the nonlinear eigenvalue problem with A = 1. 

In the other limiting case A ^> 1 the nonlinear eigenvalue problem will not 
have solutions with A = — 1. Indeed, in this case the potential V is always deep 
with the depth > A 2 regardless of the choice of fj s , so the lowest eigenvalue 
of the operator in Eq. (3.17) will be |A| ~ A 2 3> 1 (assuming that V varies 
on the length scale of order 1). This means that the solution in the form of 
the static spike AS exists only when A < 1. Note that this result was also 
obtained in [51] by the method of topological shooting. 

The numerical solution of the nonlinear eigenvalue problem shows that for 
A = 1 and arbitrary A [recall that the solutions for all other values of A can be 
obtained using the symmetry transformation given by Eq. (3.19)] there exists 
a unique stable solution for all fj s greater than some critical value fj* (see Fig. 
4). This can be explained in the following way. For large enough values of fj s 
the potential Vo, which can always produce a localized state, dominates in the 
total potential V. As the value of fj s decreases, the effect of the potential V\ 
becomes more and more pronounced, so V gradually transforms from a single- 
well to a double-well potential. This means that with the decrease of fj s the 
function 9 will tend to localize in the minima of V\ instead of the minimum of 
Vo (see Fig. 3). On the other hand, the localization of 9 in the minima of V\ 
will in turn increase VI, since the latter is self-consistently determined by 9 [see 
Eq. (3.16)]. If one constructs the solution of the nonlinear eigenvalue problem 
iteratively, for small enough values of fj s one will find that at each step of the 
iterations the potential V is such that at the next step the distribution of 9 will 
become localized further and further away from the origin. On the other hand, 
it is easy to show that there are no solution of the nonlinear eigenvalue problem 
in the form of a pair of spikes some distance L > 1 apart. Suppose that we 
have a solution in the form of two spikes centered at £ = ±L/2, with L 1. 
Let us multiply Eq. (3.17) by 0% and integrate it over positive £. As a result, 
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Fig. 5. The values of 77(0) (a) and 8(0) (b) as functions of A for the one-dimensional 
AS obtained from the numerical solution of Eqs. (3.12) and (3.13). 

using the fact that for L 3> 1 with exponential accuracy 0(0) = 0, we obtain 
the relation A 2 J °° 9 3 fj^ = 0. This relation, however, cannot be satisfied 
since the function fj^ is positive definite for all positive £, so the solution of 
the assumed form does not exist. So, for fj s < fj* the iterative procedure will 
not converge and at fj s = fj* the solution of the nonlinear eigenvalue problem 
abruptly disappears. 

From Fig. 4 one can see that A is a monotonically decreasing function of 
f] s , with its maximum attained at fj s = fj* ~ 2.40 for A = 1 (see Fig. 4). 
According to Eq. (3.19), we have |A max | = |A(f/*)| = const/A 2 . Therefore, at 
some A = A d ~ 1 we will have A max = —1, so that for A > A d there will be 
no solutions corresponding to the one-dimensional static spike AS. Thus, at 
A = Ad there is a bifurcation of the static spike AS which results in the local 
breakdown and leads to its splitting and self-replication (see Sec. 5.1 and 6). 

The numerical solution of Eqs. (3.12) and (3.13) together with Eq. (3.15) 
confirm our conclusions about the behavior of the sharp distributions as the 
value of A is varied. Figure 5 shows the dependences of the values of 9 and fj 
at £ = on A obtained from the numerical solution of Eqs. (3.12) and (3.13). 
From this figure one can see that the solution indeed disappears at A = A d 
with the value of A d found to be A d = 1.35. Figure 6 shows the distributions 
of 9 and fj in the spike obtained from the numerical solution of Eqs. (3.12) and 
(3.13) for a particular value of A. It also shows the entire solution obtained 
by matching the sharp and the smooth distributions for the particular values 
of A and e. Notice that the distributions given by Eqs. (3.21) give a very good 
approximation to the actual solution whenever A is not in the immediate 
vicinity of A d (for example, at A < 0.8 these distributions give the solutions 
with the accuracy better than 10%). 

3.1.3 Casee 1 ' 2 < A < 1 

Let us now consider the intermediate case e 1 / 2 C A <C 1. In this case the 
results of Sec. 3.1.1 and 3.1.2 both predict for 9 and 77 in the spike 
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Fig. 6. Sharp distributions (solid line) and 7y(£) (dashed line) in a 

one-dimensional static AS, obtained from the numerical solution of Eqs. (3.12) and 
(3.13) with A = 1; (b) The full solution for 6 (solid line) and r\ (dashed line) at 
A = 1 and e = 0.01. 



£(*) = - cosh" 2 (£), Vs = % (3.22) 

and the correction to rj to be given by Eq. (3.7) with 4e 2 ^r? s /3 = e [see Eqs. 
(3.5) for A > e 1 / 2 ]. Note that Eqs. (3.22) were also obtained in [51]. 

Observe that the procedure presented in Sec. 3.1.2 is valid with accuracy e 
when A ~ 1. For these values of A it is justified to assume in the match- 
ing condition that with this accuracy 77(0) = 0. As the value of A decreases, 
the actual value of 77(0) ~ e/A 2 grows, so the accuracy of the above men- 
tioned approximation decreases, and at A ~ e 1 / 2 ~ Af, this approximation 
becomes invalid. On the other hand, in the procedure discussed in Sec. 3.1.1 
the matching condition uses the true value of 77(0) ~ 1 for A ~ A b) but neglects 
the variation of 77 in the spike. The latter gives the corrections of order A 2 to 
the solution given by Eq. (3.3), which are of order e when A ~ Af, and grow as 
A increases. For A ~ e 1 / 4 both these procedures give the same solutions with 
the accuracy e 1 / 2 , so in fact for e <C 1 one can construct the solution in the 
form of the static spike AS asymptotically for all values of A. When A < e 1 / 4 , 
one should use the procedure described in Sec. 3.1.1 and when A > e 1 / 4 one 
should use that of Sec. 3.1.2. 

Above we presented the ways to construct asymptotically the solution in the 
form of the one-dimensional static spike AS. As such, these procedures should 
be good only for sufficiently small values of e <C 1. According to the analysis 
above, this AS exists in a wide range A b < A < A4 with Af, ~ e 1 / 2 1 and 
Ad ~ 1. This implies that in order for the whole asymptotic procedure to be 
in quantitative agreement with the actual solution we must have Af, <^ Ad- In 
view of Eq. (3.6), this will be the case when e < 0.03. 
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3.2 Three-dimensional radially- symmetric static spike autosoliton 



Let us now study the higher dimensional static spike ASs. These are radially 
symmetric spikes of large amplitude and the size of order e [9-11]. As we will 
show below, in the Gray-Scott model the properties of the solutions in the 
form of the higher-dimensional radially-symmetric spike ASs turn out to be 
different from those of the static spike AS in one dimension. 

Let us consider a three-dimensional AS first. The distributions of 9 and r\ in 
the form of the AS will be determined by Eqs. (2.15) and (2.16) in which 
the time derivatives are put to zero and only the radially symmetric part of 
the Laplacian is retained, with the spike centered at zero. When e < 1, one 
can once again use singular perturbation theory and separate the sharp dis- 
tributions (inner solutions) in the spike from the smooth distributions (outer 
solutions) away from the spike. 

As in the case of the one-dimensional AS, away from the spike the activator 
and the inhibitor become decoupled, so that 9 = there and the smooth 
distribution of rj is given by 



V(r) = 1 - — , (3.23) 
r 

where r is the radial coordinate and a is a certain constant [see Eq. (2.16)]. 
The constant a is determined by the strength of the 5-function like source 
term at r = 0. Integrating Eqs. (2.15) and (2.16) over the spike region, we 
obtain that 



a = - j r 2 9(r)dr, (3.24) 
o 

where the integration was extended to the whole space since 9 = away from 
the spike. One can see an important difference between Eq. (3.23) and Eq. 
(3.10) for the one-dimensional AS: in the case of the three-dimensional AS the 
derivative of i] with respect to r in the smooth distribution becomes singular 
at r = 0. This means that the scaling of 9 and r] in the spike will be different 
from that of the one-dimensional AS. Indeed, for r ~ e the value of r\ must 
be positive, so we must have a ~ e. According to Eq. (3.24), this implies that 
in the spike 9 ~ Ae~ 2 . Also, from Eq. (3.23) one can see that near the spike 
i] varies by values of order 1 on the length scale of e. Since we must have 
A9 2 i] ~ 9 in Eq. (2.15) in the spike region to have a solution, the scaling for 
the variables and the parameter A will be the following 
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0max ~e ?7min ~ 1, A ~ 6. (3.25) 

Introducing the scaled quantities and the stretched variable 



= e6, 77 = 77, i = e _1 A, f = -, (3.26) 



e 



and retaining only the leading terms in Eqs. (2.15) and (2.16), we can write the 
equations describing the sharp distributions (inner solutions) of the activator 
and the inhibitor in the spike as 



d 2 9 d-ld§ 7 ~ 2 _ ~ n , nnrrS 

ie + —dt + A( "i- e ^ (3 - 27) 



d fj d — 1 di] 

de + ~Tdt 



6 2 f) = 0, (3.28) 



with d — 3 [for generality we put an arbitrary dimensionality of space d in 
Eqs. (3.27) and (3.28)]. The boundary conditions are neutral at £ = 0, and 
zero for 9 at infinity. The precise boundary condition for fj at infinity which 
ensures the proper matching between the sharp and the smooth distributions 
has to be specified. To do this, we note that, according to Eq. (3.23), to order 
e we have 77 = 1 at r ^> e (or £ ^> 1). This means that the boundary condition 
for 77 in Eq. (3.28) must be taken to be 77(00) = 1. 

It is convenient to perform the following change of variables 



= y, V = V. + 1 j ? ; (3-29) 

where we define 7/(0) = %(0) = 0. Obviously, f] s must satisfy < fj s < 1. In 
these variables, we can write Eq. (3.27) in the form of the nonlinear eigenvalue 
problem 



d 



2 



di 



+ v(0 



9 = X9, (3.30) 



m = (n, + 2Z£) , (3.3!) 

and Eq. (3.28) as 



% = |A|0. (3.32) 
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The problem now has a one-dimensional form similar to the one in Sec. 3.1, 
but is defined for £ > 0, with zero boundary condition for 6 at £ = 0. As in 
the one-dimensional case, the solution that corresponds to the AS must be 
the lowest bound state and have A = — 1. The latter is achieved by adjusting 
the value of fj s . Also, according to the definition of 77, the matching condition 
77(00) = 1 corresponds to the boundary condition 77^(00) = 1 — fj s for Eq. 
(3.32). Integrating Eq. (3.32) with A = — 1 over £, we transform it into an 
integral condition 



00 

J Qd£ = 1 - fj s . (3.33) 





This condition fixes the normalization of 9. 



It is possible to show that the nonlinear eigenvalue problem possesses a con- 
tinuous symmetry generated by 



db 
d\ 

db 

dfj a 
db_ 

dB_ 

db 
dfj 

db 

dA 
~db 



2A, 

277,(1-77,), (3.34) 
= £(1-277,), 

?7(1 - 277 S ), 
-.-1(1-2%). 



From these equations one can see that if there is a solution of the nonlinear 
eigenvalue problem with certain f] s , A, and A, there is also a solution with 



A' = A^% A' = A 



n l/2 



(3.35) 



with 77^ arbitrary. Since X'(fj' s ) is a monotone function of 77^ that goes from 
to infinity as ^ changes from to 1, for any 77, it is always possible to choose 
a unique value of 77, for which A' = —1. So, as in the one- dimensional case, 
there is a one-to-one correspondence between the solutions of the nonlinear 
eigenvalue problem with arbitrary A and the sharp distributions. 



The potential V can be written as a sum of two parts: V — V + V\, where 
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Fig. 7. Dependence X(fj s ) for the three-dimensional radially-symmetric AS. Results 
of the numerical solution of the nonlinear eigenvalue problem for different values of 
A. 



V 



A 2 9r] s 



A 2 9{f) - 9) 
l 2 ' 



(3.36) 



The qualitative form of these potentials for sufficiently small values of A co- 
incides with the one shown in Fig. 3 for £ > 0. According to Eqs. (3.32) and 
(3.33), at these values of A we have the following estimates for V and V\ when 
fj s is close to either 1 or 



V ~-fj a (l-fj,), Vi (1 - 77,) 2 . 



(3.37) 



In writing these estimates, we used the fact that the characteristic length scale 
of the variation of 9 and fj is 1. 



3.2.1 Case A ~ e: autosoliton collapse 

For A ~ 1, one can analyze the solutions of the nonlinear eigenvalue problem in 
the following way. First of all, for sufficiently small values of A the potential V 
will be so shallow that there will be no bound states in the eigenvalue problem 
at all. When the value of A is increased, at some A = A the potential V will 
become capable to localize a state at fj s ~ \ [see Eq. (3.37)]. When the value 
of A is further increased, the minimum value of A = A min will decrease, so that 
at some A = Ab we will have A m i n = — 1 (see Fig. 7 in which the numerical 
solution of the nonlinear eigenvalue problem for several values of A is shown) . 
For A > Ab we will have A m j n < —1, so there will be two values of fj s at 
which A = — 1 (Fig. 7). Therefore, the solutions of the nonlinear eigenvalue 
problem with these values of fj s will correspond to the sharp distributions we 
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are looking for. Thus, the solution in the form of the three-dimensional static 
spike AS exists only for A > A^. From the numerical solution of Eqs. (3.27) 
and (3.28) we obtain that in this case A b = 5.8. 

It is not difficult to see that the solution corresponding to the largest value 
of fj s will be unstable. Indeed, if the value of fj s is decreased, we will have an 
increase in V\ and a decrease in Vq. Since V\ can localize 9 easier than Vq, a 
small decrease of fj s will produce such a deviation of 9 that will [through Eqs. 
(3.32) and (3.36)] further distort the potential V in the same manner. In other 
words, if we construct an iterative map that takes V, calculates the solution 9 
of the eigenvalue problem, and then generates the new V by solving for f/, it 
will take us from the unstable solution with greater fj s to the stable solution 
with lower fj s if fj s is decreased at the start, or to the trivial solution 9 = 
and fj — 1, which is obviously stable, if the value of fj s is increased. Thus, the 
solution corresponding to the stable radially-symmetric static AS should be 
unique. 



3.2.2 Case A > e: annulus 

The numerical solution of Eqs. (3.27) and (3.28) shows that for A not far 
from Aj, the distribution of 9 in the AS has the form of a spike centered at 
zero. As the value of A increases, the shape of the AS changes. To see how 
the AS behaves as the value of A is increased, a special treatment of the 
case A ^> 1 is needed. When A becomes large, the potential V contains a 
large factor of A. This factor can be compensated by choosing, for example, 
1 — fj s ~ A~ 2 <C 1, which will correspond to the unstable solution with larger 
fj s . Alternatively, one could have the potential V shifted along £, so that it is 
centered around £ = R 1. In that case the main contribution to V will be 
given by Vq ~ —A 2 /R for fj s not close to either or 1. Since 9 exponentially 
decays at large distances from £ = R, for R 1, the boundary conditions 
at £ = become inessential and can be moved to minus infinity. In this case, 
if one neglects the terms of order 1/R in the potential, one can solve the 
nonlinear eigenvalue problem exactly. The solution will be given by 



£=^^cosh- 2 (^jA, (3.38) 



4 

with 



R= A2 ^ 1 - f i'\ (3.39) 
6 



So far, we obtained a continuous family of solutions parameterized by i] s in 
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the case A >> 1. This result may seem surprising, since for A ~ 1 we showed 
that there should be only one stable solution to the nonlinear eigenvalue prob- 
lem. However, as we will show below, all the solutions found in the preceding 
paragraph except for a single solution are in fact structurally unstable, so the 
stable solution is indeed unique for A 3> 1. 

The reason for the structural instability of the solutions with R 3> 1 is that in 
the limit R — > oo the problem possesses translational symmetry. As a result, 
for R 3> 1 there is a degenerate mode corresponding to the translations of 
the spike as a whole along £. One should therefore study the stability of the 
solution with respect to that mode. 

To analyze the stability of the solutions with R ^> 1, we need to calculate the 
1/R correction to the solution obtained above. Let us write Eqs. (3.27) and 
(3.28) in the form that is valid to order 1/R 



% + c6s + + A9 2 fi -9 = 0, (3.40) 
R 

%-^ 2 = 0, (3.41) 

where for the solution sought we must have c = 0. We wrote Eq. (3.40) so 
that it reminisces an equation describing the solution traveling with constant 
speed c. 

We can use the expression for 9 obtained above to calculate the variation of 
fj in the spike. Integrating Eq. (3.41) with 9 from Eq. (3.38) and using the 
boundary condition fj^(— oo) = 0, we obtain 



71 = 71 ^2R~ l 2mcosh ^— 

+itanh 2 ^ + (3.42) 

Substituting this expression for rj into Eq. (3.40), multiplying it by 9%, and 
integrating over £, for fj s < 1/2 we obtain the following expression for c 



1 

°~ ~A 2 



2A* 



+ 



12 



i 12R f\ 24R 

V 



(3.43) 



In writing this equation we assumed the relationship between i] s and R from 
Eq. (3.39). 
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Fig. 8. The values of 77(0) (a) and 6(0) (b) as functions of A for the three-dimensional 
radially-symmetric AS obtained from the numerical solution of Eqs. (3.27) and 
(3.28). 

The analysis of Eq. (3.43) shows that we have c = only for R = R* (and 
fj s = fj*), where 

** = |, « = |. (3.44) 

Thus, the corrections of order 1/R destroy the solutions with R^ R* . 

Consider the flow generated by the equation dR/dt = c(R). The behavior of 
this flow near the fixed point R = R* determines the stability of the solution. 
According to Eq. (3.43), we have c > for R < R* and c < for R > R*, so 
the flow is into the fixed point. Therefore, the solution with R = R* is stable. 
One can also write an equation similar to Eq. (3.43) in the case fj s > 1/2, 
which corresponds to another branch of the solutions obtained above. The 
analysis of this equation shows that for those solutions c < for all fj s , so the 
flow transforms the solution into the trivial solution 9 = 0. Thus, the solution 
that corresponds to the radially-symmetric static AS is indeed unique even 
for A > 1. 



3.2.3 Comparison of the two cases 

From the arguments given above it is clear that when the value of A is increased 
from A b ~ 1 to A >> 1, the solution in the form of the AS should gradually 
transform from the spike to the annulus of large radius R ~ A 2 . This is what 
we see from the numerical solution of Eqs. (3.27) and (3.28). Figure 8 shows 
the dependence 7/(0) and 6(0) as a function of A in the radially-symmetric 
AS obtained from this solution. From this figure one can see that 7/(0) indeed 
approaches fj* = 1/3 as I increases. The dependence of the radius of the 
annulus was also found to be in good agreement with Eq. (3.44) for large 
values of A. 

Figure 9 shows the distributions of 6 and fj in the radially-symmetric three- 
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Fig. 9. (a) Sharp distributions (solid line) and t/(£) (dashed line) in a 

three-dimensional radially-symmetric static AS, obtained from the numerical so- 
lution of Eqs. (3.27) and (3.28) with A = 10; (b) Construction of the full solution 
for 9 (solid line) and rj (dashed line) at A = 0.1 and e = 0.01. 

dimensional AS for a particular value of A which is intermediate between the 
spike and the annulus. This figure also shows the entire solution in the form of 
an AS obtained by matching the sharp and the smooth distributions for the 
particular values of A and e. 

So far we studied the situations when R ^> 1 but still smaller than the in- 
hibitor length, which in these units is of order e -1 . When R reaches this value, 
Eqs. (3.27) and (3.28) will no longer be justified for the description of the 
distributions of 6 and 77 in the annulus. In the unsealed variables this will 
happen when A ~ e 1 / 2 ~ A^ [see Eq. (3.44)], where A'jp is the minimum 
value of A at which the one-dimensional static AS exists (see Sec. 3.1.1). At 
these values of A the radially-symmetric AS can be effectively considered as 
a one-dimensional AS, so when A reaches some value A d ~ e 1//2 , the solution 
in the form of an annulus will transform into a quasi one-dimensional AS of 
infinite radius. Note, however, that this bifurcation point is essentially differ- 
ent from the bifurcation at A = Ad of the one-dimensional AS (Sec. 3.1.2) in 
that it does not involve local breakdown or splitting of the AS. 



3.3 Two-dimensional static spike autosoliton 

Let us now turn to the two-dimensional case. The analysis of the two-dimensional 
radially-symmetric static AS turns out to be analogous to that of the three- 
dimensional AS, so we will not go into much detail here, but will only give the 
main results. 

As in the case of the three-dimensional AS, away from the spike 9 = and 
the distribution of r\ is given by 

ri(r) = l-aK (r), (3.45) 
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where Kq is the modified Bessel function, r is the radial coordinate, and a 
is a certain constant. The value of a is determined by integrating Eqs. (2.15) 
and (2.16) with the time derivatives set to zero over the spike region. In two 
dimensions a is given by 



a = 1 J r9(r)dr. (3.46) 
According to Eq. (3.45), when r becomes of order e, we have 

r] ~ 1 - a In e _1 + a ln(r/e). (3.47) 

In order for 77 to remain positive, we must have a ~ l/lne _1 . Then, on the 
length scales of e the variation of 77 in the spike will be of order a <C 1. In 
the same way as in the case of the three-dimensional AS, this results in the 
following scaling for the main parameters of the AS 



#max~e\ ?7min ~ T^T' A^elne 1 . (3.48) 



In the scaled variables 



6 = e9, 77 = 77 In e" 1 , A = — ^ = -, (3.49) 

erne -1 e 

the equations for the sharp distributions in the spike will take the form of 
Eqs. (3.27) and (3.28) with d = 2. Since the variation of r) in the spike is of 
order 1/ lne -1 , according to Eq. (3.47) we must have a = 1/ hie -1 in the limit 
e — > 0. This gives us the condition for matching the sharp and the smooth 
distributions. Let us introduce the variables 



= A- 1 9, fj = f) + A- 1 6-f ]s , (3.50) 

where fj s is a constant that is chosen so that 77" (0) = 0. Then, the matching 
condition can be written in the integral form as 



00 

je^=\X\-\ (3.51) 



where A is a constant that must be equal to -1. Also, in terms of the new 
variables the equation for the sharp distribution of fj can be written as 
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|A|0 (3.52) 



with %(0) = 0. 



As with the three-dimensional AS, in two dimensions the problem of finding 
the sharp distributions can be written as a nonlinear eigenvalue problem 



--M4) +m ^ x ^ (3 - 53) 

V(O = -A 2 9(Vs + V-0), (3.54) 
with the potential V that can be separated as V — V + VI, where 



V = -A 2 6r] s , V l = -A 2 6{ri - 9). (3.55) 

These potentials have the form shown in Fig. 3 for £ > 0. The lowest bound 
state with A = — 1 will give us the solution we are looking for. This condition 
is achieved by adjusting the value of fj s . The nonlinear eigenvalue problem is 
invariant with respect to the transformation given by Eq. (3.19). 

When A < 1, the potential V acquires a small factor (as in the case of 
the one-dimensional AS), which can be compensated only by choosing fj s ~ 
AT 2 3> 1. Therefore, the potential V will be dominated by V which can 
always localize a bound state with A = — 1. Notice, however, that in order for 
the approximations made to derive the equations for the sharp distributions 
to remain valid, we must have fj s < lne -1 , so in fact this argument is valid 
only down to A ~ (lne" 1 ) -1 / 2 . It is easy to show that, similarly to the one- 
and three-dimensional cases, the solution in the form of the two-dimensional 
AS will disappear at A < A b ~ e(lne _1 ) 1/2 . 

At A sufficiently small, the AS looks like a spike with the maximum value of 
9 centered at £ = 0. As in the case of the three-dimensional AS, when the 
value of A is increased, the AS gradually transforms into an annulus of radius 
R, which grows with A. According to Eq. (3.51), when R increases, we have 
9 ~ R^ 1 1 and fj ~ 9, so the potential Vo ~ A 2 /R starts to dominate. The 
analysis similar to that for the three-dimensional AS shows that for A ^> 1 
the parameters of the AS are given by 



These results are also supported by the numerical solution of the equations for 
the sharp distributions. The dependences of fj(0) and 9(0) on A are presented 
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Fig. 10. The values of fj(0) (a) and 9(0) (b) as functions of A for the two-dimensional 
radially-symmetric AS obtained from the numerical solution of Eqs. (3.27) and 
(3.28). 




Fig. 11. Sharp distributions (a) and fj(£) (b) in a two-dimensional radi- 

ally-symmetric static AS, obtained from the numerical solution of Eqs. (3.27) and 
(3.28) with A = 5. 

in Fig. 10. The solution of Eqs. (3.27) and (3.28) in the form of the two- 
dimensional radially-symmetric static spike AS at a particular value of A is 
also presented in Fig. 11. Of course, when R ~ e _1 , the approximations made 
in the derivation of the equations for the sharp distributions are no longer 
valid, so, as in the case of the three-dimensional radially-symmetric AS, the 
two-dimensional radially-symmetric AS of radius R* will transform into a quasi 
one-dimensional AS of infinite radius. According to Eq. (3.56), this will happen 
when A > A d ~ e 1/2 me -1 . 

Finally, we note that the small parameter of the singular perturbation ex- 
pansion in the two-dimensional case turned out to be 1/lne -1 , so one should 
expect it to give a good quantitative agreement with the actual solutions only 
for extremely small values of e. Nevertheless, the leading scaling given by Eq. 
(3.48) (up to the logarithmic terms) should be in good agreement even for 
not very small e. Also, it is not difficult to modify the theory in such a way 
that it uses e as a small parameter in the expansion. In that case the sharp 
distributions will contain a weak logarithmic dependence on e. 
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4 Traveling spike autosolitons 



Up to now, we only considered the solutions of Eqs. (2.8) and (2.9) that 
correspond to the static ASs. At the same time, when a = tq/t v <C 1 and e is 
sufficiently large, there exist solutions that propagate with a constant speed 
without decay — the traveling ASs [9-11], so now we are going to look for 
these solutions. Throughout this section we assume that a C 1. 

As we will show below, in the Gray-Scott model the traveling ASs are realized 
for sufficiently small a and have the shape of narrow spikes of high amplitude 
which strongly depends on a. To analyze the traveling spike ASs, it is conve- 
nient to measure length and time in the units of I and r e , respectively. Then, 
the equations describing the AS traveling with the constant speed c along the 
x-axis will take the form 



g + * + ^_, = , (4.1) 

r 2 pl + q-'c^ - e''ri + 1 - , = 0, (4.2) 

dz z az 

where we introduced a self-similar variable z = x — ct. The solution with 
c > travels from left to right. The distributions of 9 and r] should go to the 
homogeneous state 9 h = and r\ h = 1 [Eq. (2.11)] for z — > ±oo. 



4-1 Non- diffusive inhibitor: e a 



1/2 



There are two qualitatively different types of traveling spike ASs in the Gray- 
Scott model. First we consider the ultrafast traveling spike AS, which is real- 
ized when the inhibitor does not diffuse, that is, when L = (or e = oo). Such 
an AS was recently discovered by us in a similar reaction-diffusion model (the 
Brusselator) [44]. A remarkable property of this AS is that it has the shape of 
a narrow spike whose velocity c ~ AcT 1 ! 2 is much higher than the character- 
istic speed 1/tq (which in these units is of order 1) determined by the physical 
parameters of the problem, and whose amplitude goes to infinity as a — > 0. 



4-1.1 Case a 1 / 2 Ci4Ca 1 ^ 2 : ultrafast traveling spike autosoliton 

If we assume that 6 ^> 1, we can drop the last term from Eq. (4.1) and neglect 
the last two terms in Eq. (4.2) (with the term involving the second derivative 
of i] dropped in the limit of large e) in the front of the spike where 77 ~ 77^ = 1. 
If we then multiply the latter equation by A and add it to the former equation, 
we will get 
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This equation can be straightforwardly integrated. If we introduce the vari- 
ables 



9 = aA~ 1 9, c = a 1/2 A- 1 c, £ = A a - 1/2 z, (4.4) 
we can write the solution for r\ as 



~ ld9 , N 

(4-5) 

where we took into account the boundary condition 9(+oo) = 0^(+oo) = 0, 
i](+oo) = 1. Substituting this expression back to Eq. (4.1) (with the last term 
dropped), we arrive at the following equation 



,H: iW f ~c - -J 2 ) + 9 2 - 9* = 0. (4.6) 



de v c 

One can see that in Eq. (4.6) all the a- and A-dependence is absent, so Eq. 
(4.4) (with all the tilde quantities of order 1) in fact determines the scaling of 
the main parameters of the traveling spike AS for a <C 1. As was expected, the 
AS will have the speed which diverges as a — > 0. Also, note that the width of 
the front of the AS, which is of order a 1 / 2 A -1 goes to zero as a — > 0. Thus, the 
distributions of 9 and i] in the front of the ultrafast traveling spike AS will be 
given by the "supersharp" distributions (in the sense that their characteristic 
length scale is much smaller than 1) described by Eqs. (4.5) and (4.6). 

Let us take a closer look at Eq. (4.6). This equation has the form of an equation 
of motion for a particle with the coordinate 9 and time £ in the potential 
U — ^- — with the nonlinear friction with the coefficient c— 4# 2 . Since the 

3 4' c 

derivative of the friction coefficient is positive for all c, the friction increases 
as c grows, so there are no special features associated with its nonlinearity. 
For 9 > the potential U has a maximum at 9 — 1 and a minimum at an 
inflection point 9 = 0. The supersharp distribution of 9 will therefore be the 
heteroclinic trajectory going from 9 = 1 to 9 = 0. 

It is clear that if the friction is not strong enough, the particle starting from 
9 = 1 will miss the point 9 = and go to minus infinity, so we must have 
c > c*, where c* is some positive constant of order 1. On the other hand, it is 
clear that when c > c*, the particle will always get from 9 = 1 to 9 = 0, so in 
fact there is a continuous family of such solutions. Thus, we have a multiplicity 
of the front solutions and, therefore, a selection problem [5]. To answer the 
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question about the front selection, we need to consider higher-order corrections 
to the solution of Eq. (4.6) coming from Eqs. (4.1) and (4.2). According to 
these equations, for small 9 the next order correction will amount to adding 
the term — aA~ 2 9 to Eq. (4.6). In this situation the potential U will actually 
have a maximum at 9 = and a minimum at 9 m i n ~ aA~ 2 , so only the 
trajectory with the minimum velocity c will reach 9 = 0, whereas all other 
trajectories will be stuck at 9 = 9 m i n . Thus, we can conclude that in the limit 
a — > the selected front solution in our problem has the velocity c = c*. The 
numerical solution of Eq. (4.6) shows that the value of c* is c* = 0.86. The 
numerical simulations of Eqs. (2.8) and (2.9) confirm these conclusions. The 
main parameters of the traveling spike AS, therefore, are 



#max = Aa-\ c = 0.86 x Aa- 1/2 . (4.7) 

Note that the numerical solution of Eq. (4.6) in the form of the supersharp 
front differs from ss h = |[1 — tanh(0.50£)] by less than 1%. Also note that 
the results given by Eq. (4.7) precisely coincide with those obtained by us for 
the Brusselator [44] . This is due to the fact that the supersharp distributions 
in these two models are described by the same equations. 

In the back of the supersharp front the value of 9 goes exponentially to 1, and 
rj goes exponentially to [see Eq. (4.5)]. Note, however, that in writing the 
equations describing the supersharp distributions we neglected the last two 
terms in Eq. (4.2). When the value of 77 decreases, at rj ~ a 2 A" 2 the term 9 2 r\ 
becomes of order 1, and the equations for the supersharp distributions cease 
to be valid. This will happen at a distance of order o^^A^ 1 In AcT 1 behind 
the location of the supersharp front. We can therefore call the region of this 
size right after the front where r\ exponentially decays to some value 77 min the 
secondary region of the supersharp distributions. Since the width of this region 
is still much smaller than 1, we can assume that 9 = # max there. Then, the 
distribution of r/ in the secondary region of the supersharp distributions is 
given by Eq. (4.2) in which we should drop the last term, since i] <^ 1 there. 
We obtain 



Vssh2 = a 2 A~ 2 + Ce^, (4.8) 

where the constant C should be determined by matching to the asymptotics 
of the supersharp distribution of i] at £ — > —00 (this requires an explicit 
knowledge of the solution in the supersharp region). As can be seen from this 
equation, we have ?7 m i n = a 2 A~ 2 . 

As z passes the secondary region of the supersharp distributions, 9 2 r] becomes 
of order 1, and therefore can be dropped from Eq. (4.1). Then the activator 
and the inhibitor become decoupled, so the characteristic length scale of the 
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variation of 9 significantly increases. According to Eq. (4.1), for c ^> 1 the 
characteristic length scale of the decay of 9 behind the supersharp front is 
of order c ~ Aa^ 1 ^ 2 , which is still much smaller than the length scale of the 
variation of rj behind the spike (the refractory region), which is of order oT x c 
(see below). This means that after the secondary region of the supersharp 
distributions we should find the region of the sharp distributions. According 
to Eq. (4.1) with the terms d 2 6/dz 2 and A9 2 i] dropped, the solution for 9 in 
this region will be 



9 sh {z) = a- l Ae z ' c , (4.9) 

where we chose the position of the supersharp front to be at z = (with the 
accuracy of a). This expression for 9 sh can be substituted back into Eq. (4.2) 
to calculate i]^. The analysis of this equation then shows that one can neglect 
both a~ 1 c d7]/dz and —r] in the region of the sharp distributions, so i]^ and 
# s h are related locally. The resulting expression for the sharp distribution of rj 
takes the following form 



Vshl = a 2 A- 2 e~ 2z / c . (4.10) 

As will be shown in the next paragraph, this equation is in fact valid only in 
the part of the sharp distributions region, so we will call it the primary sharp 
distribution of rj. 

According to Eq. (4.10), the value of 77 exponentially grows behind the region of 
the sharp distributions, so at some distance of order a~ 1 / 2 A\na~ 1 A 2 one can 
no longer neglect the term a~ x c dr]/dz in Eq. (4.2). If we take this derivative 
into account, we can solve Eq. (4.2), provided that 6 is still given by Eq. (4.9). 
The solution will have the following form 



W*) = |r(0,ef(— o))e e§( ^ o) , z = ^\n2aA- 2 , (4.11) 

where F(a,x) is the incomplete gamma function. In writing the last equation 
we matched this solution with the one from Eq. (4.10) at large z — z . We will 
call this distribution of rj the secondary sharp distribution. 

For yet more negative values of z the distribution of i] approaches rj ~ — ac^z 
[see Eq. (4.11)], so the characteristic length scale of the variation of rj becomes 
of order cT l c 3> c. This means that we enter the refractory tail of the AS 
where r\ relaxes to rjh, that is, the region of the smooth distributions. For 
these values of z the distribution of 9 already relaxed to zero, so Eq. (4.2) can 
be easily solved. To do that we should recall that up to c or x c the region 
of the sharp distributions is located at z — 0, and to the leading order in a 
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Fig. 12. Distributions of 9 and r] in the ultrafast traveling spike AS: (a) the back 
of the AS; (b) the front of the AS. Results of the numerical solution of Eqs. (2.15) 
and (2.16) with e = oo, a = 1CT 3 , and A = 1. Length is measured in the units of I. 

we have 77(0) = 0. This immediately gives us the solution for 77 in the region 
of the smooth distributions 

Vsm = l-e az / c . (4.12) 

The entire solution in the form of the ultrafast traveling spike AS is presented 
in Fig. 12. This figure actually shows the result of the numerical simulations 
of Eqs. (2.15) and (2.16) with a = 1CT 3 and A = 1. One can see an excellent 
agreement of this solution with the distributions obtained above. 

Thus, we introduced an asymptotic procedure for constructing the solution 
in the form of the traveling spike AS in the Gray-Scott model in the limit 
a — > 0. This solution is considerably different from the solutions in the form 
of the traveling ASs in N-systems (see Sec. 1). In N-systems the speed and 
the distribution of 9 in the AS front are determined only by the equation for 
the activator with 77 = r\ h in the limit a — > 0, so the speed of the AS cannot 
exceed the value of order 1 [9,11]. The distribution of 9 in such an AS can 
be separated into two regions: the region of the sharp distributions, which 
corresponds to the moving domain wall whose characteristic size is of order 
1, and the region of the smooth distributions, where the distribution of 9 is 
slaved by the distribution of 77 which varies on the length scale of oT x [11,27]. 
In other words, in the limit a — > there is only one boundary layer in the 
solution for 9 (within a single domain wall) and no singularities in the solution 
for 77. In contrast, in the Gray-Scott model the speed and the amplitude of the 
ultrafast traveling spike AS become singular in the limit a — > 0. Moreover, 
there are three regions with different behaviors for 9 in the ultrafast traveling 
AS: the region of the supersharp distributions, where 9 varies on the length 
scale of a 1 / 2 , the region of the sharp distributions, in which the characteristic 
length scale of 9 is a -1 / 2 , and the region of the smooth distributions, where 
= 0. The latter happens to be a specific property of the Gray-Scott model, 
in more general models the distribution of 9 is slaved by the distribution of 77 
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Fig. 13. (a) Solution of Eqs. (4.14) and (4.15) for A = 3.76. (b) Dependence c(i) 
obtained from Eqs. (4.14) and (4.15). 

in the smooth distributions and thus has the characteristic length scale of its 
variation of order a~ l c [44]. Moreover, the distribution of i] can be separated 
into five regions where the asymptotic behavior of rj is different. In other 
words, the solution in the form of the ultrafast traveling spike AS contains 
four boundary layers in the limit a — > 0. 

4-1-2 Case A ~ a 1 / 2 : disappearance of solution 

According to the procedure presented above, the main parameters of the AS, 
such as the amplitude and the velocity are determined solely by the supersharp 
distributions of 9 and 77. However, according to Eq. (4.7), when A becomes 
of order a 1 / 2 , the velocity of the AS becomes of order 1, so the separation of 
the distributions of 9 and 77 into the supersharp and the sharp distributions 
in the spike becomes invalid. For these values of A the treatment of the spike 
region has to be modified. Note that according to Eq. (4.7) we still have 
#max ~ or 1 ! 2 ^> 1 for A ~ a 1 / 2 . Let us introduce the following variables 



9 = a 1/2 9, fj = rj, A = a - l/2 A. (4.13) 
In these variables we can write Eqs. (4.1) and (4.2) as 

9 zz + c9 z + A9 2 f)-9 = 0, (4.14) 

cfj z -9 2 f) = 0, (4.15) 



where we neglected the last two terms in Eq. (4.2). These equations with the 
boundary conditions 6*(±oo) = and f/(+oo) = 1 can be solved numerically. 
Figure 13(a) shows the solution of these equations for a particular value of A. 
One can see that the distribution of 9 has the form of an asymmetric spike, 
while the distribution of rj goes from fj — 1 at plus infinity to fj — i] m i n at 
minus infinity. The numerical solution of Eqs. (4.14) and (4.15) shows that 
the traveling solution exists only for A > A^t = 3.76. The numerical solution 
also shows that for A > A bT we have r) min < 0.05, which decreases with the 
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increase of A, so with good accuracy we can assume that r] min = 0. Figure 
13(b) shows the dependence c(A) obtained from the numerical solution of 
Eqs. (4.14) and (4.15). Observe that already for A > lAA^ the velocity of 
the AS agrees with Eq. (4.7) with accuracy better than 15%. Behind the spike 
region (in the region of the smooth distributions) we will have 9 = and 
r\ — 1 — (1 — 77 m i n ) e az / c . If one assumes that i] min = 0, one would arrive once 
again at Eq. (4.12). 

We would like to emphasize that even for A ~ a 1 / 2 <C 1, that is, for such 
a small deviation of the system from thermal equilibrium, the amplitude of 
the AS is 6* max ~ or 1 ! 2 ^> 1 and the velocity c ~ 1. In other words, the AS 
remains a highly nonequilibrium object in the system only slightly away from 
thermal equilibrium. 



4-1.3 Case A ~ a 1 ^ 2 : oscillatory tail 

In the other limiting case A ^> 1 the behavior of the secondary sharp distri- 
butions in the back of the AS acquires special features. This is due to the fact 
that for A ^> 1 the phase trajectory 9{rj) in the phase plain of 9 and 77 may 
pass close to the unstable fixed point [see Eq. (2.12)] 9h3 — A, r\ya — A~ 2 , so 
the behavior of the distributions of 9 and 77 behind the spike becomes oscilla- 
tory, so 9 and rj will not be able to get back to the homogeneous state 9 = 9^ 
and 77 = rjh at z — —00. To see that, let us introduce 

9 = a ^ 2 9, rj = a^rj, A = a 1/2 A, £ = -. (4.16) 

c 

Then, we can rewrite Eqs. (4.1) and (4.2) behind the spike as follows 



9 6 + A9 2 f] - 9 = 0, (4.17) 
% - 9 2 fj + 1 = 0, (4.18) 

where we kept only the leading terms. In order for the solutions of these 
equations to properly match with the primary sharp distributions, we must 
have 9 ~ e^ and fj ~ e~ 2 ^ as £ — > +00 [see Eqs. (4.9) and (4.10)]. The numerical 
solution of Eqs. (4.17) and (4.18) with these boundary conditions shows that 
at A > A dT = 0.90 the distributions of 9 and fj become oscillatory behind the 
spike. Thus, we conclude that the ultrafast traveling spike AS exists in a wide 
range A bT < A < A dT , where A bT = 3.76a 1/2 and A dT = 0.90a~ 1/2 . Notice 
that the oscillatory behavior of the distributions of 9 and 7] behind the spike of 
the AS is essentially related to the Hopf bifurcation of the homogeneous state 
= 9h3, V = Vh3 for 0.41a -1 / 2 < A < oT x l 2 (see Sec. 2). On the other hand, for 
A > a -1 / 2 this homogeneous state becomes stable, so in that case the traveling 
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spike AS transforms to a wave of switching from one stable homogeneous state 
to the other. 

4-1-4 Justification of the condition e 3> a 1 / 2 and case a = e 2 

Above we considered the case, in which the inhibitor does not diffuse. Let us 
see how the diffusion of the inhibitor affects the properties of the ultrafast 
traveling spike AS. Since the main parameters of the AS are determined by 
the primary supersharp distributions, the diffusion of the inhibitor should 
not significantly affect these distributions. According to Eq. (4.2), the term 
e~ 2 d 2 n/dz 2 ~ e~ 2 a~ 1 A 2 is small compared to the leading terms which are 
of order a~ 2 A 2 [see Eq. (4.4)] if e 2 3> a regardless of AP] In terms of the 
physical parameters of the problem this means that the ultrafast traveling AS 
will be described by the solution obtained above as long as D e ^> D v , where 
Dq and D v are the diffusion coefficients of 9 and n, respectively. It is also clear 
that when a ~ e 2 , the properties of the AS will not change qualitatively. An 
interesting special case a = e 2 which corresponds to the activator and the 
inhibitor with equal diffusion coefficients can be treated in an analogous way 
(see also [52]). The resulting equation for the supersharp distributions will have 
the form of Eq. (4.6), but without the nonlinear friction term. This equation 
can be solved exactly, giving in this case the velocity c = l/y/2 [52,53], which 
is in fact close to the one obtained in the case of the non-diffusing inhibitor. 
The explicit expression for the supersharp distributions in this case are: 9 SB h = 
\ (l -tanh^), Vsshl = § (l+tanh^), and n ssh2 = a 2 A~ 2 + e^A The 
rest of the solution will be the same as in the case e = oo. All this allows us to 
conclude that the ultrafast traveling spike AS in the Gray-Scott model exists 
when e > a 1 / 2 . 

4-2 Diffusive inhibitor: e <C a 1 / 2 

Now let us analyze the second type of the traveling spike AS which is realized 
when both a <C 1 and e <C 1. It is obvious that in this situation there exist 
a solution in the form of the spike AS whose velocity is equal to zero (see 
Sec. 3.1). What we will show below is that when a becomes small enough, in 
addition to the static spike AS there are solutions in the form of the traveling 
spike AS which propagates with constant velocity whose magnitude is > l/r e . 

Since e <C 1, it is natural to separate the distributions of 9 and n into the sharp 
and the smooth distributions. As in the case of the one-dimensional static AS, 
in the spike region we introduce the scaled variables from Eq. (3.11), with 
£ = z in this case. In terms of these variables, Eqs. (4.1) and (4.2) become 

2 Note that these estimates remain valid even at A ~ Ay?. 
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Fig. 14. Qualitative form of the traveling spike AS at e 2 <C a < e for c ~ 1 (a) and 
c > 1 (b). 



22 + c# 2 + A6 2 f) -6 = 0, 
fizz - 2 fj = 0, 



(4.19) 
(4.20) 



where we assumed that e <C a 1 / 2 and 9 2r q ^> 1 in the spike and kept only 
the leading terms. Similarly to the one-dimensional case, the scaling for the 
variables in the spike region is given by Eq. (3.8). According to Eq. (4.20), the 
asymptotic behavior of fj at large \z\ is 



fj ~ K±Z, Z — > ±00, 



(4.21) 



where k± are some constants. Therefore, the distributions of 9 and fj in the 
spike will qualitatively have the form shown in Fig. 14. 

It is convenient to introduce the variables from Eq. (3.14). Then, Eq. (4.20) 
can be written as 



Vz 



9 - c9 7 . 



(4.22) 



Because of the translational invariance, we have the freedom to choose the 
position of the spike. We will do it in such a way that the maximum of 9 is 
located at z — 0, that is, we have 9 Z (0) = 0. Also, according to Eq. (4.22), 
we can add an arbitrary function a + bz to its solution, so we may replace 
fj — > fj + fj s + kz, where fj s and k are arbitrary constants, and require that 
77(0) = 77^(0) = 0, so the function fj(z) is completely determined by 9(z). In 
view of all this, Eq. (4.19) becomes 



9 ZZ + c9 z + A 2 9 2 (f] s + Kz + rj-9)-9 = 0. 



(4.23) 



According to Eq. (4.21), we have f/ 2 (±oo) = k±—k. This implies that fj z (+oo) — 
fj z (— 00) = k+ — K-. Integrating Eq. (4.22) over z, we obtain an integral rep- 
resentation of this condition 
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+oo 

J 9dz = k + — k_ 

— oo 



(4.24) 



By a similar integration, the value of k is determined as 
o 

k = K-+ J 9dz. (4.25) 

— oo 

To study the solutions of Eqs. (4.22) and (4.23), we use the mechanical analogy. 
Equation (4.23) can be considered as an equation of motion for a particle of 
unit mass with the coordinate 9 and time z in the potential U = — y + 
A 2 fj S Y — A 2 ^ in the presence of friction, with the friction coefficient c, and an 
external time-dependent force — A 2 9 2 [kz + fj(z)]. For fj s > 2A~ l the potential 
has two maxima at 9 = and 9 = 9 m , and a minimum in between. The 
particle slides down from the maximum of the potential at 9 = and after an 
excursion toward 9 = 9 max < 9 m at z = returns to 9 = 0. Notice that since 
by definition fj(0) — fj z (0) — 0, the value of fj should be rather small in the 
spike region, so one could think of the second term in — A 2 9 2 [kz + fj(z)] as a 
small perturbation. 

In the presence of friction the external time-dependent force acting in Eq. 
(4.23) must be such that it accelerates the particle at some portion of the 
trajectory. According to Eq. (4.22) and the fact that 9 Z < for z > 0, we have 
fj > for those values of z. Since the values of k relevant to our analysis are 
positive (see below), in the portion of the trajectory where z > the external 
force does accelerate the particle. All the kinetic energy that is gained by 
the particle during this part of the motion must be dissipated by the friction 
force, so that the particle arrives at 9 = with zero velocity. This defines 
the precise value of the friction coefficient c as a function of fj s and k. Recall 
that in addition the distribution 9(z) must satisfy the integral condition in 
Eq. (4.24), so in fact we are not free to choose the value of fj s . Thus, for given 
values of k± there may exist only a discrete set of the velocities c. 

Away from the spike 9 is zero, and r] is described by the smooth distributions. 
If we introduce the variable ( = ez, we can write Eq. (4.2) in the form 

?7 CC + P' 1 ^ + 1 - 77 = 0, (4.26) 

where (3 = a/e. We should choose such a solution of this equation that cor- 
rectly matches with the sharp distribution. To do that, we use the fact that 
to order e the value of 77(0) = 0, so the smooth distribution of rj is 
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V(C) = 1 -exp(-/c±C), 



(4.27) 



where 



k± = c±v|+W (4 28) 

and one should take k + if £ > 0, or k_ otherwise. Note that for /3 ~ 1 and 
c ~ 1 we have k± ~ 1. From Eq. (4.27) one can see that when £ approaches 
the spike region, we have r\ ~ re ± £ = eK±z. This means that we should use the 
values of k± given by Eq. (4.28) in solving the inner problem given by Eqs. 
(4.22) - (4.25). 



4-2.1 Case eA 2 < a <C eA: bifurcation of the static and the traveling au- 
tosolitons 

Equations (4.23) - (4.25) are difficult to deal with and in general can only be 
treated numerically. Such a treatment was recently performed by Reynolds, 
Ponce-Dawson, and Pearson, who also derived these equations in a different 
context [50]. The problem can be significantly simplified in the case A< 1. 
In this case there is a small factor of A 2 multiplying a number of terms in 
Eq. (4.23). It can be partially compensated by choosing fj s ~ A" 2 ^> 1. If we 
neglect the other terms containing A 2 and put c = in Eq. (4.23), we can 
solve this equation [together with Eq. (4.24)] exactly. The result is 



to^^cosh- »(£), = (4.29) 

Note that according to Eqs. (4.25) and (4.29), we may put k = (k + + k_)/2 = 
c/(2[3). 

The equation describing the small deviations 89 = 9 — 6q due to the order A 2 
terms that were dropped in the derivation of Eq. (4.29) can be obtained by 
the linearization of Eq. (4.23) around 9 . Assuming that c< 1 and retaining 
only the term that gives the nontrivial contribution to c, we get 



59 = A 2 9 2 kz + c-^. (4.30) 

dz 

The operator in the left-hand side of this equation has an eigenvalue zero, 
corresponding to the translational mode 59 = d9 /dz. Therefore, in order for 
Eq. (4.30) to have a solution, its right-hand side must be orthogonal to this 
translational mode. The operator in Eq. (4.30) is self-adjoint, so in order to get 



dz 2 



2A%9 + 1 
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the solvability condition for this equation, we should require that the integral 
of its right-hand side multiplied by dd^/dz be equal to zero. With the use of 
Eq. (4.29), this gives us the velocity c as a function of k± 



A 2 \ 



K 2 , — K 2 



C = 



(4.31) 



To determine the velocities that are actually realized in the traveling spike 
AS, we need to take into account the matching conditions that are given by 
Eqs. (4.28). With the use of these equations, Eq. (4.31) becomes 

Since in the derivation of this equation we assumed that c 1, it will be valid 
only for (3 <C A. 

Two observations can be made from Eq. (4.32). First, at some (3 = (3 T (A) [or 
at some A = At(/3)] the velocity of the AS becomes zero. Since this happens 
for (3 ~ A 2 and a ^> e 2 , we have e 1 / 2 < i C 1, so there is also a solution 
with c = (see Sec. 3.1). The presence of a point where the velocity of the 
traveling solution goes to zero therefore signifies a bifurcation of the static 
AS. Second, according to Eq. (4.32) the velocity of the obtained solution is a 
decreasing function of A (or an increasing function of f3). In contrast, we would 
expect the velocity of the traveling spike AS to be an increasing function of 
the excitation level A. 

Let us consider an iterative map that takes c, substitutes it to Eq. (4.28) with 
the fixed k + + K-, calculates k and uses Eq. (4.31) to give the new value of c. 
Clearly, those c given by Eq. (4.32) (or c = 0) are the fixed points of this map. 
However, an analysis of this map shows that an arbitrarily small deviation of 
c from that given by Eq. (4.32) will lead to the unlimited growth of c if the 
deviation is positive, or to zero if the deviation is negative. In other words, the 
fixed point given by Eq. (4.32) is unstable. Also, one can easily see that for 
A < At the fixed point at c = is stable for A < At (or (3 > f3r) and unstable 
otherwise. This means that the solution with non-zero velocity we found above 
and the static solution at A > A T or (3 < f3 T should be unstable. Therefore, 
the stable traveling solutions should have the velocity c > 1 and may exist 
both when A < At and A > At- Also, when A > At, the solutions with c = 
should be unstable, so the static spike AS spontaneously transforms into the 
traveling spike AS, whose speed c > 1. These conclusions are also supported 
by the numerical simulations of Eqs. (2.8) and (2.9). 
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4-2.2 Case a <C eA: ultrafast traveling autosoliton 

Above we considered the situation in which c C 1. Let us now study another 
possibility: c > 1. In this case the distribution of 9 will become strongly 
asymmetric [see Fig. 14(b)]. Indeed, according to Eq. (4.23), we will have 
~ or cz at z — > +oo and 9 ~ e z / c at z — > — oo. In other words, we can once 
again separate the distributions of 9 and fj into the regions of the supersharp 
distributions (in the front of the spike) and the sharp distributions (in the back 
of the spike). In the region of the supersharp distributions the supersharp front 
will have the width of order c" 1 <C 1. Let us introduce £ = cz. Then, we can 
write Eq. (4.22) integrated over £ and Eq. (4.23) as 



9tt + 9t + c- 2 A 2 9 2 (rj s + rj-9) = 0, (4.33) 
% = -9 + C, (4.34) 

where we retained only the leading terms and moved the point where 9 at- 
tains its maximum value to minus infinity [this amounts to putting n = 
in Eq. (4.33) and redefining the boundary condition for fj to be fj(— oo) = 
fj^(-oo) = 0]. One can see that by rescaling 9 and fj with # max all the A- and 
^max-dependence from Eqs. (4.33) and (4.34) can be absorbed into c. These 
equations have a solution in the form of a supersharp front connecting 9 = 
ahead of the front with 9 = 9 max behind the front only when fj s = # max . In this 
case we have %(+oo) = 9 max [see Eq. (4.34)], what corresponds to K- — and 
k+ = cfi~ x . According to Eq. (4.28), these values of k± can only be realized 
when c ^> /3, with 9 max = f3~ l . Note that the fact that fj s = # max behind the 
supersharp front means that fj exponentially decays to zero at £ — > — oo. 

The numerical solution of Eqs. (4.33) and (4.34) shows that the velocity of 
the supersharp front is 



c = 1.22 x A/T 1 . (4.35) 

Since by assumption c ^> 1, we must have /3 <C A. Recall that in the derivation 
we also assumed that c» (3. According to Eq. (4.35), this leads to (3 2 A. 
Since, as we will show in Sec. 4.2.3, the solution in the form of the traveling 
AS exists only when /3 < 1, this condition is always satisfied when f3 A. 
Note that the numerical solution of Eqs. (4.33) and (4.34) in the form of the 
supersharp front differs from ^ ss h = |[1 — tanh(0.44£)] by less than 0.5%. 

In the region of the sharp distributions the characteristic length scale of the 
variation of 9 is c, so we can neglect the term 9 ZZ in Eq. (4.23). Recalling that 
fj = in this region, we can write the solution of this equation that represents 
the sharp distribution of 9 behind the supersharp front as ^ s h = 9 ma , x e Z//c , 
where we assumed that the supersharp front is located at z — 0. 
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Fig. 15. Qualitative form of the dependence c(A) (a) and the dependence c{f3) (b) 
for the traveling spike AS. 

When the value of A is decreased, the velocity of the unstable traveling solution 
grows and the velocity of the stable traveling solution decreases until they 
reach the value of order 1 when the approximations used in deriving the above 
equations become invalid. According to Eq. (4.35), this will happen at A ~ (3, 
so at some A = A^t ~ (3 the solution in the form of the traveling spike AS 
will disappear. Therefore, the velocity of the traveling spike AS as a function 
of A or (3 should have the form shown in Fig. 15. 

For f3 < 1 the traveling AS exists at A > f3. On the other hand, for these 
values of (3 the static spike AS remains stable up to the values of A ~ ft 1 / 2 
[see Eq. (4.32) and the discussion below]. Therefore, for (3 < A < (3 1 ^ 2 the 
static spike AS will coexist with traveling. 



4-2.3 Case A ~ e 1 , a < e: oscillatory tail 

When the speed of the traveling spike AS increases, the behavior of the dis- 
tributions of 9 and r\ in the back of the AS changes in the way similar to the 
case of the ultrafast traveling spike AS (Sec. 4.1.3). For large enough values of 
A the sharp distributions in the back of the spike become oscillatory. To see 
that let us introduce the variables 



e = e 1 ' 2 9, rj = e-\ A = e 1 ' 2 A, c = e^c, £ = -, (4.36) 



where c is given by Eq. (4.35). Then, keeping only the leading terms, we can 
write Eqs. (4.1) and (4.2) in the region of the sharp distributions as follows 



% + A6 2 f) -6 = 0, (4.37) 

+ - ~° 2f i " 1 = °- ( 4 - 38 ) 

To match the solutions of these equations with the supersharp distributions 
we must have rj — > for £ — > +oo, so 9 ~ e^ and f\ ~ e~ 2 ^ for large £. 
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Fig. 16. Dependence (3d(A) for the traveling spike AS. 

Similarly, in order for the distributions of 9 and rj to properly match with the 
smooth distributions, for £ — > — oo we must have 0(— oo) = and 77^ = — /3 
[see Eq. (4.28) with c >> 1]. Equations (4.37) and (4.38) with these boundary 
conditions can then be solved numerically. As a result, we find the values of 
(3 = (3 d as a function of A at which the behavior of the distributions of 9 and 
fj in the back of the spike changes to oscillatory. The dependence (3d(A) is 
plotted in Fig. 16. From this figure one can also see that the solution in the 
form of the traveling spike AS exists only at (3 < (3 C ~ 1.6. The analysis of 
Eqs. (4.37) and (4.38) also shows that for A 3> f3 3 ^ 2 the second derivative of 
fj in Eq. (4.38) can be dropped, so the problem reduces to that given by Eqs. 
(4.17) and (4.18) with A replaced by Af3 1 ^ 2 . Therefore, for large enough values 
of A we should have (3d — 0.81A -2 . The numerical simulations of Eqs. (2.15) 
and (2.16) show that the onset of the oscillatory behavior in the back of the 
spike in the traveling spike AS results in the splitting of the AS as it moves 
(see Sec. 6). 

4-2.4 Case a ~ e, A ~ 1: qualitative considerations 

So far, we concentrated on the situations in which either c < 1 or c > 1. We 
found that only the solutions with c ^> 1 should be stable and can be realized 
only if f3 < 1 and <C A (Sec. 4.2.2). As was shown in Sec. 4.2.3, in the case 
A ^> 1 the solution in the form of the traveling spike AS will exist only if 
P < P c ~ 1. It is clear that qualitatively these conclusions can also be made 
when both A ~ 1 and f3 ~ 1, what corresponds to c ~ I. In this case the 
AS will not be able to propagate at A < A^t ~ 1 and will transform into the 
static AS. On the other hand, for A > Adr ~ 1 the traveling spike AS will 
start splitting. When the value of (3 is increased, the value of A^t will grow 
and the value of A d r will decrease, so at some (3 — (3 C the solution in the form 
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of the traveling spike AS will disappear. This conclusion is supported by the 
numerical simulations of Eqs. (2.15) and (2.16), with the value of f3 c found to 
be very close to the one obtained in the preceding paragraph (Sec. 6). 

The reasons for stopping of the traveling spike AS at A < A bT or splitting 
at A > A dT are the following. Consider Eqs. (4.22) and (4.23), using the 
mechanical analogy. Let us see what happens as the value of A decreases 
when c ~ 1. In this case the external force — A 2 9(kz + fj) accelerating the 
particle on the way from 9 max to becomes smaller [see Eq. (4.23)], while the 
friction force remains the same. When the value of A becomes small enough, 
the dissipation will exceed the acceleration, so the particle will not be able to 
reach 9 = and the solution traveling with constant velocity will disappear. 
On the other hand, as long as A b < A < A d , where A b is the point where the 
solution in the form of the static one-dimensional AS disappears and A d the 
point where the static one-dimensional AS starts splitting (see Sec. 3.1.2), the 
static solution will exist, so the traveling AS will be able to stop and transform 
into static when the value of A is decreased. 

According to the definition of 9, when z is close to zero, the time-dependent 
external force is dominated by — A 2 9 2 kz, which is positive in some interval 
zq < z < 0. This means that the particle is accelerated by this force before it 
stops at 9 — 9 max . If the value of A is big enough and the friction coefficient 
c ~ 1, the forces from the potential U and the friction may not be enough to 
stop the particle before it reaches the maximum of the potential at 9 — 9 m . 
Then, if the particle moves past 9 m , it will keep on moving toward plus infinity 
and will never be able to return to 9 = 0. This means that the solution in the 
form of the traveling spike AS must disappear at some A = A d r ~ 1- Notice 
that the same argument can be applied to the static AS at A = A d , so this 
bifurcation point should indeed correspond to the onset of splitting. 



5 Stability of static spike autosolitons 

In Sec. 3 we showed that in the Gray-Scott model the solutions in the form 
of the static spike ASs exist in certain ranges of the parameter A in arbitrary 
dimensions. According to the general qualitative theory of ASs [9-11], these 
static ASs are stable in a wide range of the system's parameters. However, 
when these parameters are varied, the ASs may become unstable with respect 
to various kinds of instabilities [9-11]. So, now we are going to study the 
stability of these solutions against different types of perturbations in the full 
system given by Eqs. (2.15) and (2.16), with e<C 1. 
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5.1 Static one- dimensional autosoliton 



We start by analyzing the one- dimensional static spike AS in higher dimen- 
sions. The equations describing small deviations 59 = 9 — 9 and 5i] = 77 — 770 of 
the activator and the inhibitor, respectively, from the distributions Oq(x) and 
770 (x) in the form of the static one-dimensional AS are obtained by linearizing 
Eqs. (2.15) and (2.16) (we chose the axes so that 9 and i] depend only on 
x). Let us take 



SO = 59 kul (x) e 



iut—iky 



Sr] = Sr) kuj (x)e 



iujt—iky 



(5.1) 



where u is the complex frequency and k is the wave vector that characterizes 
the transverse perturbations of the AS. Then, after some algebra, Eqs. (2.15) 
and (2.16) linearized around 9 and rj can be written as 



d 2 



dx 2 



+ 1 + iuj + k 2 - 2A9 r] 



d 2 

-e~ 2 — + 1 + ia 1 tu + t' 2 k 2 
ax 1 



59 ku) = A9 2 5r] ku) , 
Srj ku ,= 



(5.2) 



-A' 1 



d 2 



dx 2 



+ 1 + iu + k z 



59 kLU , 



(5.3) 



where we measured length and time in the units of I and Tg, respectively. 
Equation (5.3) can be solved by means of the Green's function 



Srjku, = 



-59 



ku) 



e (1 + iu — ie 2 a 1 uj) 
2AVI + m^u + e~ 2 k 2 



+00 



^ lx - x ' 1 S0Ux')dx', 
(5.4) 



where we neglected a term of order e 2 . This expression for 5i] kw can be sub- 
stituted back into Eq. (5.2) to get a single operator equation in terms of 59 ku) 
alone. Note that the operator in the left-hand side of Eq. (5.2) is a Schrodinger 
type operator with the potential in the form of the well of the depth and the 
size of order 1. Indeed, according to the results of Sec. 3.1, the characteristic 
length scale of the variation of 9q is of order 1, and we have #0 ~ Ae^ 1 and 
rjo ~ A~ 2 e in the spike, so A9o1]q ~ 1 and exponentially decays away from the 
spike. Also, the function 5r] kuJ in the right-hand side of Eq. (5.3) is multiplied 
by the function 9\ which also exponentially decays away from the spike, so in 
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fact we are only interested in the behavior of S6k w in this region. Therefore, 
in solving the operator equation mentioned above one should only know the 
sharp distributions of 9 and r\ in the spike. If we take these distributions writ- 
ten in terms of the variables from Eq. (3.14), we can reduce Eqs. (5.2) and 
(5.3) to 



d 2 



dx'' 



+ l + iu + k — 2A%9 - 2A 2 9 fj + 3A Z 9< 

+00 



59. 



e-M¥ 2 (l + icu - ie 2 q-y) 
2y/l + ia~ l uj + e~ 2 k 2 



c 



ey/l+i a -iu+e-W\x-x'\ S9 kul ( X ')dx' . 

(5.5) 



This is the basic equation for studying the stability of the static one- dimensional 
AS in higher-dimensional systems, which is asymptotically valid for e C 1. One 
has to solve this equation as an eigenvalue problem, find the modes 89 n and 
the values of u = to n {k) corresponding to them. The instability of the AS will 
occur when the real part of 7 = —iuj is negative.0 To analyze the stability of 
the static spike AS in one dimension, one should simply put k = in Eq. (5.5). 
Note that since the right-hand side of Eq. (5.5) is an exponentially decaying 
function of x, a mode 59 n can be unstable only if it is localized since otherwise 
we would have 7 > 1 + k 2 > 0. In other words, the unstable modes should be 
in the discrete spectrum of the solutions of Eq. (5.5). 



5.1.1 Dangerous modes 

To solve Eq. (5.5) in general is a formidable task. Since we are mostly inter- 
ested in the instabilities, we can simplify the problem by identifying the most 
dangerous modes and try to solve Eq. (5.5) for these modes. In order to get an 
idea as to what these dangerous modes might be, we first consider Eq. (5.5) 
with k = and uo = 0. In this case Eq. (5.5) can be written asymptotically as 
an eigenvalue problem 



d 2 



dx 2 



+ 1 - 2A%0 O - 2A%fj + 3A 2 9 2 



69 r , 



e~ x A 2 9l 



+00 



59 n (x')dx = X n 59 n , (5.6) 



where the unstable modes 59 n will correspond to X n = 0. This equation deter- 
mines the instabilities of the spike AS in one dimension for sufficiently large 



3 Re 7 is the damping decrement (decay rate) of a fluctuation. We will use u and 
7 interchangeably throughout this section. 
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values of a. Note that these modes (in the PDE sense) should be in correspon- 
dence with the bifurcation points of the solution O and 770 (in the ODE sense). 
Also note that equations of the type of Eq. (5.6) were analyzed by Kerner and 
Osipov in the context of the stability of the spike AS in systems of small size 
[9-11,16]. 

When Aj, < A <C 1, one can neglect the last two terms inside the brackets 
in Eq. (5.6). Since the solution O and fj s in this region is known [see Eqs. 
(3.3) and (3.5)], the potential in the Schrodinger operator becomes V(x) = 
1 — 3 cosh~ 2 (:r/2). The eigenfunctions of this operator are well-known [54]. The 
spectrum of this operator contains three discrete eigenvalues 



Ai 0) =--, * = J- cosh" 3 *, 
4' V 32 2' 

\P = 0, 50? ] = Mtmh J cosh- 2 f , (5.7) 
V 8 2 2 



x (o) _ 3 



S e® = ^5tanh 2 I - l) cosh" 1 1, 



and a continuous spectrum of the eigenvalues A^ = 1 + k 2 > 1, character- 
ized by a wave vector k. The shape of the functions from Eq. (5.7) is shown 
in Fig. 17. The integral operator in Eq. (5.6) has a discrete eigenvalue A s = 

|2 

with the eigenfunction 50 s = | cosh- 4 (x/2) oc 0% with 



1 + J1 A2 



A 2 

^ r 1 v " A2 » 

the adjoint 89* = 1, and a degenerate continuous spectrum of the eigenfunc- 
tions satisfying <50(x)cfc = (and for the adjoint ^ 2 (x)^*(a;)dx = 0) 
which all correspond to A = 0. 

If the value of A were much smaller than Aj, (ignoring for the moment the fact 
that 9q is defined only for A > Aj, and thinking about it as A- independent), 
the integral operator in Eq. (5.6) would be small, so the solution of Eq. (5.6) 
would be determined by the the spectrum of the Schrodinger operator [Eq. 
(5.7)]. When the value of A is increased, at A ~ A b the integral operator starts 
mixing these eigenmodes. From the form of the eigenfunctions 59^ and 59 s it 
is clear that most of the mixing will occur between 59^ and 59 s . It is also clear 
that because of the above mentioned properties of this integral operator Eq. 
(5.6) [as well as Eq. (5.5)] should generally have a discrete spectrum containing 
three different eigenvalues with the eigenfunctions that look like 59^\ 59f\ 
and <5#2°\ respectively. Therefore, we will denote these functions as 59q, 59i, 
and 502, respectively. The fluctuation 50 o will result in the variation in the AS 
amplitude, the fluctuation 50\ in the shift of the AS along the x-axis, and 502 
in the widening of the AS. 

It is not difficult to see that Eq. (5.6) with the above mentioned simplifications 
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Fig. 17. Fluctuations 66^' (a), 59^' (b), and 59^> (c) from Eq. (5.7). 

is satisfied for 59$ = cosh~ 2 (:r/2) and Ao = at A = A^. This is in agreement 
with the fact that the solution in the form of the static spike AS disappears 
at A = A b . Also, this equation has a solution with Ai = and 56i = 56^. 
This is in fact the consequence of the translational symmetry of the problem 
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[9-11]. In general, Eq. (5.6) is always identically satisfied with Ai = for 
50i = d9~o/dx, which for A ~ A b is up to a coefficient given by 56^. Note 
that this fact alone does not correspond to the actual instability, for which 
we should have Re 7 < 0, so a careful analysis of the full Eq. (5.5) is needed 
to determine whether the function S9\ is in fact a dangerous mode. Also note 
that by symmetry 89\ is completely decoupled from S9o and 59 2 . 

When A b <C A <C 1, we will have |A S | ~ AA 2 /Al 3> 1, so in this case it is the 
Schrodinger operator that becomes a perturbation to the integral operator in 
Eq. (5.6). Therefore, the mode 59 = 59 s with A ^ A s can no longer give an 
instability. On the other hand, we still have the mode 59 2 that for A ^> A b 
should be orthogonal to 89 s . Since 59 2 is not strongly coupled with 89 s , it is 
natural to expect that the value of A2 should be positive for < A < 1 
(recall that for A = in Eq. (5.6) we would have it equal to A 2 °^ = 3/4 > 0), 
so this mode should not give an instability for these values of A either. 

In contrast, when A becomes of order 1, one should also consider the last two 
terms in the integral operator in Eq. (5.6). These terms will give a negative 
contribution to A2 (see Sec. 5.1.2), so at some A ~ 1 it will become zero, 
signifying an instability of the static spike AS in one dimension with respect 
to the 59 2 mode. Comparing this with the results of Sec. 3.1, we conclude that 
this will happen precisely at A = A d = 1.35 at which the solution in the form 
of the static one-dimensional spike AS disappears. Thus, we would expect that 
the modes £#0,1,2 are the only dangerous modes in the whole range of existence 
of the static spike AS. 



5.1.2 Case a ^> 1 and k = 0: stability of the autosoliton in one dimension 

Let us now go back to the analysis of Eq. (5.5). In the case a»l and k = 
it can be asymptotically written as 



When A ~ A b , Eq. (5.8) can be further simplified by noting that in this case 
its right-hand side is of order 1 and the last two terms in the square bracket in 
its left-hand side can be neglected. Using 9 and fj s from Eqs. (3.3) and (3.5), 
we can write Eq. (5.8) in the case A ~ A b asymptotically as 



dx 2 



+ 1 - 7 - 2A%9 - 2A 2 9 f] + 3A 2 9 2 59 




(5.8) 
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Fig. 18. Dependence Re 70,2 C<4) for the static one-dimensional AS for A ~ (a) 
obtained from the numerical solution of Eq. (5.9), and A ~ Ad (b) obtained from 
Eq. (5.8). 
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cosh- 4 (I) J 56{x')dx'. (5.9) 



Naturally, Eq. (5.8) in the case i < 1 should give the same results as Eq. 
(5.9) in the case A > A b . 

To see that the AS is actually stable in one dimension for Ab < A < Ad in the 
case of sufficiently large a we solved Eqs. (5.8) and (5.9) numerically, using 
the asymptotic distributions 9 and fj obtained numerically in Sec. 3.1.2. 
We solve these equations by discretizing the operators and diagonalizing the 
obtained matrices (taking the limit e — > in the case A ~ 1). This gives 7„ 
as functions of A. The numerical solution indeed shows that for A b < A < A d 
the damping decrement Re 7„ is positive for all 56 n and that it goes to zero 
for 59q as A — > A b or for S62 as A — > Ad signifying the instabilities of the AS 
at these points. The shape of the fluctuation 56$ suggests that the instability 
at A = A b will result in the collapse of the AS, while the shape of 89 2 suggests 
that the instability at A = A d will result in the local breakdown in the AS 
center and the onset of splitting. 

The damping decrements Re •jq^A) obtained from the numerical solution of 
Eqs. (5.8) and (5.9) are shown in Fig. 18. Notice that it turns out that only 
at A b < A < 1.0QA b or at 0.90 < A < A d we have Re u; ,2 = 0, so that there 
are two distinct modes 5#o,2 with different values of 7. For all other values of 
A we have Re u 7^ and there are two eigenfunctions: 59 2 and its complex- 
conjugate, which correspond to the two complex frequencies uj and — u*. For 
A b A A d we have Re 7 ~ 1 and Re u ~ 0.15. 

Observe that Eq. (5.9) can be analyzed rigorously. This analysis turns out 
to be rather involved, so we present it in appendix A. It agrees with the 
conclusions of this section concerning the case A< 1. 
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5.1.3 Case a < 1 and k = 0: instability with respect to the pulsations 

General qualitative theory of ASs suggests that for small enough values of 
a the static spike AS may become unstable with respect to the fluctuations 
with Re uj 7^ 0, resulting in the onset of the AS pulsations [9-11,43,45]. The 
analysis of Eq. (5.5) shows that in order for it to have a solution with Im uj = 
and Re uj = uj ^ we must have e 2 < a < 1 and uio ~ 1. Let us introduce 
k, 2 = e 2 /a < 1. Then, dropping 1 in the square roots in the right-hand side of 
Eq. (5.5), and putting k = 0, we get asymptotically 



In order for the exponential in the right-hand side of this equation to decay at 
large distances, we must choose the analytic branch of the square root that has 
a positive real part for all values of ujq, what is achieved by making a branch 
cut along the positive imaginary axis. Obviously, the solutions of Eq. (5.10) 
will come in pairs, for each solution with uj there will be another solution with 
—uj* . Note that because of the choice of the branch cut this equation cannot 
have solutions with Im uj > for which Re uj = 0. Therefore, an instability 
of the AS described by Eq. (5.10) must necessarily be a Hopf bifurcation (the 
exception here is the S9i mode). 

To find the instabilities of the static spike AS with respect to the pulsations, 
we solved Eq. (5.10) numerically with uj = ujq, with ujq real. This numerical 
solution of Eq. (5.10) shows that the static spike AS indeed becomes unstable 
with respect to the 59 mode with Re uj = ujo(A) at a < a w (A) ~ e 2 when 
A ~ 1. The plots of a u (A) and ujq(A) are shown in Figs. 19(a) and 19(b), 
respectively. Alternatively, for a given value of a ~ e 2 the static spike AS 
becomes unstable with respect to the pulsations when A < A u < 1. Note that 
the onset of the instability of the static spike AS with respect to the pulsations 
in the Gierer-Meinhardt model was studied by Osipov by the multifunctional 
variational method [43]. 

Equation (5.10) can be simplified in the case At, <ti A 1 and e 2 < a <C 1. 
In this case one can neglect the last two terms in the left-hand side of Eq. 
(5.10), the exponential and the term ik 2 ojq in its right-hand side, and use the 

4 Strictly speaking, the branch cut should begin at uj = ia [see Eq. (5.5) with 
k = 0] , which for sufficiently small a can be considered to be at zero in the analysis 
of most of the dangerous modes. This, however, cannot be done for the mode 50\ 
where the bifurcation is saddle- node (see the discussion below). 
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Fig. 19. Dependences a w {A) (a) and ujq{A) (b) for the static one-dimensional AS 
obtained from the numerical solution of Eq. (5.10); the dependences a^A) (c) and 
loq (d) obtained from the numerical solution of Eq. (5.13) at A ~ A b . 

distribution of 9 given by Eq. (3.22). As a result, Eq. (5.10) can be written 
in this case as 



^ + l + i-o-3 cosh- (| 



59 = 

A 2 (l + iu} ) cosh" 4 



x 



•/2) 



f oo 



5e(x')dx'. 



(5.11) 



— oo 



One can see that the dependence on the system's parameters enters this equa- 
tion only via the combination A 2 / k = a 1 ^ 2 A 2 /e. Solving Eq. (5.11) numeri- 
cally, we obtain that the AS destabilizes at 



a u = 5.15 x e 2 A~ 4 , u = 0.534, for A b < A < 1. (5.12) 

Note that Eq. (5.11) can be rigorously analyzed by the method similar to the 
one used for Eq. (5.9). This analysis is presented in appendix B. It yields the 
same conclusions as above and Eq. (5.12). 

The result of Eq. (5.12) is in agreement with the one obtained in [51] by a direct 
solution of Eq. (5.11). Notice that this equation gives a good approximation 
for a w [Fig. 19(a)] only for A < 0.5. Recalling that A > A b = ^/l2e [Eq. 
(3.6)], we see that Eq. (5.12) can give a good approximation only for e < 0.01 
in a limited range of A. 
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According to Eq. (5.12), as A decreases the value of increases, and when 
A reaches the value of order A b we must have ~ 1, so for these values of 
the parameters one can no longer neglect 1 compared to oT x uj in the square 
roots in Eq. (5.5). On the other hand, in this case one can still use the same 
approximations as in deriving Eqs. (5.11), except we should use Eqs. (3.3) 
and (3.5) instead of Eq. (3.22). As a result, for a ~ 1 and A ~ A b we get 
asymptotically 



dx 2 



+ 1 + iluq — 3 cosh 



3A 2 (1 



1UJ ) 



8AI VI + ia- 1 ^ 




86 = 



A 2 
A 2 



cosh 4 ( — 



56(x')dx'. 

(5.13) 



The results of the numerical solution of this equation for q u and uj are pre- 
sented in Figs. 19(c) and (d), respectively. From these figures one can see that 
for a > ao = 0.33 the considered instability of the static spike AS cannot be 
realized for any value of A. This result is in agreement with the general qual- 
itative theory of ASs [9-11]. The numerical solution of Eq. (5.10) also shows 
that all other instabilities with Re u ^ occur at significantly lower values of 
a, when the AS is already unstable with respect to the pulsations. 



5.1.4 Case a ^> 1 and k 7^ 0: instability with respect to the corrugation 

Let us now see what happens in higher dimensions when A ~ 1 and a is 
sufficiently large [so that the terms proportional to a^ 1 in Eq. (5.5) can be 
dropped] as the value of k is varied. Since the mode 861 generally requires a 
careful consideration, we will first look only at the modes 86 0t 2- If k < e, the 
operator in the right-hand side of Eq. (5.5) contains a large factor of order e _1 
and the exponential in Eq. (5.5) can be neglected. Therefore, for these values 
of k we should either have 7^ ~ 1, what corresponds to the modes at the 
bottom of the continuous spectrum (see Sec. 5.1.1), or the modes 86k should 
be orthogonal to 86 s . In the latter case the /c-dependence in the right-hand 
side of Eq. (5.5) becomes inessential, so for these k the values of 7 fc will be 
close to 7 n obtained for k = 0. Therefore, according to the results of Sec. 5.1.2, 
neither of the modes 86k should be unstable for k < t for A b < A < A d . 

The situation changes when e k < 1. Then, one can neglect 1 compared to 
e~ 2 k 2 in the square roots in Eq. (5.5) and write it asymptotically as 
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dx 2 



+ 1 - lk + k 2 - 2A%6 - 2A 2 9 fjo + 3A 2 6 2 



59 k 



2k 



+00 

J e -k\x-x'\ 56k ^ dx '_ (5. 14 ) 



Note that for k ~ 1 the coefficient multiplying the right-hand side of Eq. 
(5.14) becomes of order 1. Here we should expect a corrugation instability 
with respect to the mode 58o with some k — ko ~ 1 [9-11,45]. Indeed, as the 
value of k increases, the magnitude of the right-hand side of Eq. (5.5) decreases 
as 1/k, while the contribution to the left-hand side of Eq. (5.5) increases as 
k 2 . Therefore, for some k = k the contribution of both these two terms to 
Eq. (5.5) will be minimal, so that we can get an instability: Re 7^,, < 0. 

To show that there is indeed an instability at k ~ 1, we solved Eq. (5.14) 
numerically. Figure 20 shows the solutions for Re 7^ obtained for a particular 
value of A = 1.2. 

For k > 0.29 Eq. (5.14) has two localized solutions and a continuous spectrum 
of 7^ for A = 1.2, all having Re oo(k) = 0. The curve at the bottom of the 
figure corresponds to 59 , the curve in the middle to 59 2 , and the curve on 
the top is the bottom of the continuous spectrum. From Fig. 20 one can see 
that the AS is unstable with respect to 59q (the corrugation instability) for a 
range of wave vectors k ~ 1 . The analysis of 7& for different values of A shows 
that the static one- dimensional spike AS in higher dimensions is unstable with 
respect to corrugation for all values of A at which it exists. 

When 0.05 < k < 0.29, the complex frequency u(k) for the modes S6 2 
acquires a real part for A = 1.2. The corresponding eigenfunctions for these 
modes, as well as 7^, are complex-conjugate. For yet smaller values of k the 
real part of u vanishes once again, so we have two distinct solutions. The 
latter is related to the presence of two distinct solutions in the case of the 
one-dimensional AS studied in Sec. 5.1.2, which is obtained from Eq. (5.14) 
in the limit k — > 0. When the value of A is decreased below A = 0.90, these 
two distinct solutions disappear and the solution with a nonzero real part of 
uo(k) goes all the way to k — 0. 

Note that in the case A b < A <C 1 it can be easily shown that the AS is 
unstable with respect to the corrugation instability. Indeed, for i <C 1 and 
k ^> A 2 the operator in right-hand side of Eq. (5.14) can be considered as a 
perturbation. Therefore, in the zeroth approximation we will have 58q = 56^, 
where 69^ is given by Eq. (5.7). Then, in the first order of the perturbation 
theory we will have 
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Fig. 20. Damping decrement Re jk for the static one-dimensional AS obtained from 
the numerical solution of Eq. (5.14) for different modes at A = 1.2. 



^=(-4 + !048F + l H 1 + !048Fj • (515) 

which is negative for k < 1. Notice that according to this equation the fastest 
growing mode has k ~ 0.60, what coincides with very good accuracy with 
the results of the numerical solution of Eq. (5.14) (see Fig. 20). Similarly, for 
A b <C A 1 and a <^ it is easy to show that the AS will be unstable. 
Here we can also treat the right-hand side of Eq. (5.11) as a perturbation, 
so once again 59 = 56^. The solution of Eq. (5.11) in the first order of the 
perturbation theory with a/e 2 <^ A 4 and the exponential set to 1 then gives 
us 7 ~ — | + 27 1 Q 2 ^ A 2 a 1 / 2 e~ 1 < 0. Notice that by equating this expression to 
zero, one gets the value of q u which differs from the one in Eq. (5.12) by only 
10%. 
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5.1.5 Case a ^> 1 and k ^ 0: instability with respect to wriggling 



Let us now turn to the mode 89\. Since, because of the translational invariance 
Eq. (5.5) is identically satisfied for k = and 7 = with 59i = d9 /dx, it is of 
special interest to study its solutions for |7| <C 1 and k <C 1. The small values 
of k and 7 introduce a weak perturbation to the operators in Eq. (5.5) with 
k = and 7 = 0. Therefore, in the leading order of the perturbation theory 
we must take 59\ = d9 /dx, and in order to get 7^ we must multiply Eq. (5.5) 
by the adjoint function 59\ and integrate over x. As a result, in the first order 
in 7 fc and k 2 we obtain 



|2/i r 2„-^ + 00 + 00 



-Ik + k 2 + lk 1 ^ ' J J 9 2 (x)59* 1 (x)\x-x'\59 1 (x')dxdx' 

—00 —00 
2 +00+00 

= — (7* -cte- 2 k 2 ) [ [ 9 2 Jx)59* 1 (x)(x-x') 2 S9 1 (x , )dxdx', (5.16) 

—00 —00 

where we expanded the exponential in Eq. (5.5) up to the second order in e and 
used the normalization 59\59idx = 1. Recalling that 59i = d6 /dx, we can 
calculate the integral in the right-hand side of Eq. (5.16). Using the symmetry 
properties of 661, we find this integral to be 2 f*™ x6^56\dx f^™ 9 dx < 0, 
since for x < we have 89\ > and vice versa. By a similar integration, one can 
show that the integral in the left-hand side of Eq. (5.16) is also negative. Then, 
it is easy to see that when a > e (otherwise the AS would be unstable in one- 
dimension, see Sec. 5.1.6), for small values of k we will have 7 fc ~ — e~ 1 A 2 k 2 < 
0, so the one-dimensional static spike AS is in fact always unstable with respect 
to the 89\ mode with small k for A ^> A^. These fluctuations should lead to 
wriggling of the AS. Moreover, it is possible to show that the static spike AS 
which is stable in one dimension is in fact unstable with respect to wriggling 
for all values of A > A b . Indeed, for A ~ A b the operator in the right-hand 
side of Eq. (5.5), which is the only non self-adjoint operator there, can be 
considered as a perturbation, so we can assume that 50 1 = 591 = 59? and 
take #0 to be given by Eq. (3.3). Substituting this into Eq. (5.16), we obtain 
that for small values of k the damping decrement 7^ of the fluctuations leading 
to wriggling is given by 



7fc ^ k 2 




(5.17) 



From this equation one can see that 7^ < 0, signifying an instability, for 
A > A b and sufficiently small k. Thus, in summary, the one-dimensional static 
spike AS is always unstable in higher dimensions, so the instabilities that are 
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realized for sufficiently small a are meaningful only for the one-dimensional 
system. 



5.1.6 Case a < 1 and k = 0: instability with respect to the onset of the 
traveling motion 

In addition to the instability of the AS with respect to wriggling, another 
instability may be realized when a C 1 [45]. As was already mentioned above, 
for k = Eq. (5.16) has a simple zero solution for any a due to the translational 
invariance. However, at some special value of a = cxt the solution 7 = may 
become doubly degenerate. This will happen when the coefficient in front of 
7 in Eq. (5.16) vanishes. If this is the case, in addition to the trivial solution 
we will have another non-trivial zero solution, for which the value of 7 should 
change sign as the value of a passes through ct T . This signifies an instability 
that leads to the onset of the traveling motion. It is not difficult to see that 
this instability will occur when a < cit ~ e for A < 1. In particular, when 
A < 1, we can put 89\ = 89i = 56^, since again the only non self-adjoint 
operator is the operator in the right-hand side of Eq. (5.5) and is small. After 
a little algebra, we obtain that the instability will occur at 



Note that the last formula is in agreement with the results of Sec. 4.2.1. 
To analyze the dependence ax(A) for A ~ 1 we solved for the conjugate 
function 59\ numerically and then substituted it to Eq. (5.16). The resulting 
dependence cit{A) is presented in Fig. 21. Observe that for A < 1 the value 
of «t is given by Eq. (5.18) with an accuracy better than 10%. 



5.1.7 Comparison of the pulsation and the traveling instabilities 

According to Eqs. (5.12) and (5.18), for A > 1.58c 1 / 6 we have «t > a^, so 
if one starts with the static spike AS at A ~ 1 and sufficiently large value 
of a and then gradually decreases a, the AS will destabilize with respect to 
the fluctuation 56i and will transform into traveling. If, on the other hand, 
we have A < 1.586 1 / 6 at the start, the AS will destabilize with respect to 
the pulsations and will collapse after a few oscillations of its amplitude (see 
Sec. 6). Also, according to Eqs. (5.12) and (5.18), for a < a c = 0.83e 4//3 the 
one-dimensional static spike AS will be unstable regardless of the value of A, 
so we conclude that this AS can be excited only when a > a c . 




A b < A < 1. 



(5.18) 
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Fig. 21. Dependence ax (A) for the static one-dimensional AS obtained from the 
numerical solution of Eq. (5.16). 

5.2 Stability of the three-dimensional radially- symmetric static spike autosoli- 
ton 

So far we were studying the one-dimensional static spike AS in higher dimen- 
sions. Let us now turn to the radially-symmetric spike ASs. It turns out that 
up to the logarithmic terms, all the results on the stability of the two- and the 
three-dimensional radially-symmetric static spike ASs have the same depen- 
dence on e, so we will concentrate mostly on the three-dimensional AS. For 
the three-dimensional radially-symmetric static spike AS the small deviations 
59 = 9 — 9 and 5r] = rj — i] in the spherical coordinates can be taken as 

59 = e^Y nm ^,^)59 nu} (r), Srj = e^ Y nm ($, ^)^(r), (5.19) 

where Y nm are the spherical harmonics. Then, Eqs. (2.15) and (2.16) linearized 
around 9 and r] become 



d 2 2d n(n + l) , . n ~ aA _ ' 
dr z r dr r A 

d 2 2d n(n + 1) 9 9 



^ naJ = e- 1 ^ 2 ^, (5.20) 



dr 2 r dr 

= -eA- 1 



d 2 2d n(n + 1) 

+ „ + 1 



dr 2 r dr 



59 nuJ , (5.21) 



where, once again, length and time are measured in the units of / and r e , 
respectively, and we used the scaled variables given by Eq. (3.26). Equation 
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(5.21) can be solved by means of the Green's function, which in the case of the 
operator in the left-hand side of Eq. (5.21) multiplied by r 2 is (similar Green's 
function was used in [21,22]) 



G n0J (r,r') = < 



' Jn+i/2p?Vl+ia 1 u)K n+1/2 (erWl+ia M ^ 



rsr, 



J n+l/2(e?-Vl+ia 1 u>)K n+1/2 (tr\/l+ia ^ 



(5.22) 



r > r , 



where I„+i/2(^) and K n+1 / 2 {x) are the modified Bessel functions. Solving Eq. 
(5.21) with the use of Eq. (5.22), we obtain 



5rjnu = -eA- l 56 nw - ei _1 (l + iuo - ie 2 ^ 1 ^) ^ G nw (r, r')^(r , )r ,2 rfr'. 

o 

(5.23) 

Substituting this expression into Eq. (5.20), we obtain the following equation 



d 2 2d n(n + l) 1 . n ~~ _ ~ 2 
ar z r dr r A 



59 n 



= -6 2 {1 + ie'a-'u) J G nu) (r, r')56 nuJ (ry 2 dr' . (5.24) 

o 

This equation determines the complex frequencies of different fluctuations as 
the functions of the control parameters and has to be solved as an eigenvalue 
problem. The instability of the AS will occur when Re 7 < 0, with 7 = — \u. 



5.2.1 Case a ^> e 2 : instability with respect to the non radially- symmetric 
fluctuations 

Let us first look at Eq. (5.24) at a ^> e 2 . In this case the terms proportional 
to a~ l in the right-hand side of Eq. (5.24) can be neglected, so the Green's 
function G nul can be expanded in e. Then Eq. (5.24) becomes asymptotically 



dr 2 



2d n(n + l) , 



r dr 



+ 1 - 7„ - 2A9 fj + 



59, n 



gg(r)(l- 7 n) 
2n + 1 



g n (r,r')50 n (rydr', (5.25) 



where 
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9n{r,r 



(r/r') n , r < r' 

V ' ' ~ (5.26) 



The operator in the left-hand side of Eq. (5.25) is a Schrodinger type op- 
erator with the attractive potential — 2A9ofjo, repulsive potential d\ and the 
centrifugal potential n(n + l)/r 2 . 

As in the case of the one-dimensional AS, the modes 59 n that can lead to 
instabilities are localized. From the translational symmetry of the problem 
we know that there is a localized solution 59\ = dQ^jdr for n = 1 which 
corresponds to 7 = 0. It is clear that since the potential in the left-hand 
side of Eq. (5.25) is able to localize a state at n — 1 in the presence of 
the centrifugal potential, it will also be able to localize a state at n = in 
its absence. According to Sec. 3.2, we should in fact have an instability at 
the point A — A\> — 5.8 where the solution in the form of the static three- 
dimensional radially-symmetric AS disappears. Also, it is natural to expect 
that the operator in the left-hand side of Eq. (5.25) will also be able to have 
a bound state for sufficiently small n and sufficiently large A, when the role 
of the centrifugal potential decreases (see below). These will be the dangerous 
modes that we need to consider. 

It is not difficult to show that for A 3> 1 the AS will be unstable with respect 
to the fluctuations with n^O. Indeed, according to the results of Sec. 3.2, for 
these values of A the AS has the form of an annulus of large radius R ~ A 2 
[see Eq. (3.44)]. According to Eqs. (3.38) and (3.42), the leading contribution 
to the potential in the left-hand side of Eq. (5.25) is V(r) = —3 cosh -2 f 2 -^)- 
For R >> 1 and 1 <C n <C R all the other terms can be considered as a 
perturbation to the Schrodinger operator with this potential, so in the first 
order of the perturbation theory 50q = 89^^ — R), where 56^ is given by 
Eq. (5.7). In the first order of the perturbation theory we obtain 




This equation shows that for R 3> 1 there is a range of values of n for which 
7 n < 0. The fact that the annulus must be unstable with respect to the shape 
deformations is an obvious consequence of its quasi one-dimensional character 
(see Sec. 5.1.4). 

The presence of the localized solution 86 n of Eq. (5.25) with 7 n = for n > 
signifies an instability of the radially-symmetric static AS with respect to the 
non radially-symmetric fluctuations resulting in the distortions of the spike 
and leading to its splitting (see Sec. 7). It is clear that this instability will 
occur easier at smaller values of n, so we have to solve Eq. (5.25) for n = 2 
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(the case n — 1 corresponds to the translation of the AS as a whole without 
any deformation). The numerical solution of Eq. (5.25) shows that the AS 
indeed becomes unstable with respect to the fluctuation with n = 2 when 
A > A c2 = 8.4. In general, the numerical solution of Eq. (5.25) shows that the 
static radially-symmetric spike AS is stable only in the narrow range Ab < 
A < A C 2. It is interesting to note that Eq. (5.27) with n — 2, together with 
Eq. (3.44), gives the correct value of A c2 with accuracy better than 5%. 

When the value of A is increased beyond A c2 , the buildup of the instability 
with n = 2 will result in the splitting of the AS. After such a splitting event, the 
two newborn ASs will go apart until they are separated by a sufficiently large 
distance, and then will split again. Thus, the considered instability will result 
in self- replication of ASs (see Sec. 7). We would like to emphasize that the 
character of the instability that leads to self-replication in three-dimensional 
(and two-dimensional) systems is significantly different from that in one di- 
mension. In the former the instability results in the buildup of the shape 
distortion that eventually leads to splitting, while in the latter the instability 
leads to the widening of the activator distribution profile and local breakdown 
in the AS center (see Sec. 6). Thus, self- replication of the AS in one dimension 
is qualitatively different from that in higher dimensions. 

5.2.2 Case a < e 2 : instability with respect to the pulsations and the onset of 
the traveling motion 

According to the general qualitative theory of ASs, when a becomes small, 
the AS may become unstable with respect to the pulsations with Re uj ^ 0. 
The analysis of Eq. (5.24) with n = shows that this instability may be 
realized only when a ~ e 2 and Re to — Uq ~ 1. In view of this fact, we can 
drop 1 in the square roots in Eq. (5.22). Equation (5.24) can then be solved 
numerically. The results of this numerical solution are presented in Fig. 22. 
The upper curve in Fig. 22(a) shows the critical values of for the onset 
of this instability for those values of A at which the AS is stable at a > 6 2 . 
Figure 22(b) shows the frequency u of the fluctuations at the threshold of 
the instability. From Fig. 22(a) one can see that the static three-dimensional 
AS is unstable for all values of A if a < a c ~ 3.7e 2 . 

As was already noted, in addition to the modes studied above, we always have 
a dangerous mode 69 = d0/dr with n — 1. The analysis of Eq. (5.24) shows 
that in addition to the trivial solution, at some a = ax ~ e 2 this equation can 
have another solution with Re 7 = 0. As in the case of the one-dimensional 
AS, this solution will signify the instability which results in the AS starting to 
move as a whole. To find this instability point, we need to know the behavior 
of the Green's function Gi u (r,r') at small values of iuj. Expanding the Bessel 
functions in Eq. (5.22) and finding the adjoint function 59\ numerically, we 



63 





stable 




C : 
O : 




bo 


3 : 

- "o : 




splittii 


cfl ■ 
O : 

- a : 


pulsations 


: -7 — traveling 






Fig. 22. (a) Stability diagram for the three-dimensional radially-symmetric AS. (b) 
Frequency loq at the threshold of the pulsation instability. Results of the numerical 
solution of Eq. (5.24). In (a) the upper solid line shows the values of a^/e 2 corre- 
sponding to the onset of the pulsation instability (n = 0) , the lower solid curve 
shows the values of ax/e 2 corresponding to the instability leading to the onset of 
the traveling motion (n = 1). The dotted lines correspond to the instabilities of the 
static AS at A = A^ and A = A&. 



calculate the coefficient in front of iuj in the first order of the perturbation 
theory, keeping only the leading terms in e. The above mentioned instability 
will be realized when this coefficient vanishes. The lower solid line in Fig. 
22(a) shows the values of a-r/e 2 at which this happens. One can see that 
this instability occurs when the AS is already unstable with respect to the 
pulsations. Note that the numerical solution of Eq. (5.24) shows that the 
instabilities with respect to the fluctuations with n > 2 always happen for 
significantly lower values of a. 



5.3 Stability of the two-dimensional radially- symmetric static spike autosoli- 
ton 



Finally, we will briefly discuss the stability of the two-dimensional static spike 
AS in the limit e — > 0. In the cylindrical coordinates r, the small deviations 
56 = 9 — 9 and 5r) — rj — rjo can be taken as 



69 = S9 nu) (r) e^-^, 5 V = 5 Vnu) (r) e^"^, 



(5.28) 



where n is an integer. 



The equation for 59 nuJ is obtained by linearizing Eqs. (2.15) and (2.16) around 
9q and t)q and eliminating 5f] nuJ by inverting the equation for 5f] nuJ . As a result, 
using the scaled 9 and fj given by Eq. (3.49), we arrive at the following 
equation for 59 nu) 
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where G n (r,r') is given by 



G n (r,r') 



r < r', 
r > r', 



(5.30) 



= < 



I n {er' \/T+lcF^)K n {er\/Y+^cr^)^ 



where I n and K n are the modified Bessel functions. The analysis of this equa- 
tion can be performed in the way completely analogous to that of Eq. (5.24). 
It is clear that the instabilities of the two-dimensional static spike AS will be 
qualitatively the same as those of the three-dimensional AS. We would like 
to emphasize that as in the case of the three-dimensional AS, splitting and 
self-replication of the static spike AS in the two-dimensional system is related 
to the buildup of the fluctuation with n = 2 describing a non-symmetric dis- 
tortion of the AS. Because of the very slow convergence of the asymptotic 
theory in two dimensions (recall that the small parameter here is 1/ hie -1 ) we 
will not present a detailed study of Eq. (5.29). 



6 Pattern formation scenarios in one dimension 

In the following two Sections we present the results of the numerical simula- 
tions of the Gray-Scott model for sufficiently small e and a and compare them 
with the results of our asymptotic analysis. 

A word needs to be said about the numerical methods we used in various parts 
of the paper. Whenever it was not said otherwise, length and time were mea- 
sured in the units of L and t v , respectively. To simulate the original reaction- 
diffusion equations we used a simple explicit second-order scheme both in one 
and two dimensions. We used neutral boundary conditions in all our simula- 
tions. In order to resolve the details of the shape of the spike, sufficiently small 
spatial discretization step was needed. It was found that the step Ax = 0.25/ 
gave the solutions with accuracy of a few per cent. The stiffness of equations at 
small e or a makes the simulations very time consuming, so two-dimensional 
simulations were done on the massively parallel supercomputer (SGI-Cray Ori- 
gin 2000). A typical running time for the simulations shown in the paper is 
about 1 hour on 16 processors. 

In addition to the direct simulations of the original reaction-diffusion equa- 
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Fig. 23. (a) Distributions of 6 and r/ in a static spike one-dimensional AS; (b) 
Potential V from Eq. (3.18). Results of the numerical solution of Eqs. (2.15) and 
(2.16) with e = 0.05, a = 1, and A = 1.4. 

tions, we used various methods in solving the equations obtained from the 
asymptotic procedures. In solving for the sharp distributions in the case of 
the static spike ASs we used relaxation methods with second-order discretiza- 
tion. In solving the nonlinear eigenvalue problems we used a combination of 
the relaxation method and iterative method. In the case of the traveling spike 
ASs we used shooting method to solve the obtained ODEs. In analyzing the 
stability of the static spike ASs we used second-order discretization of the lin- 
ear operators in a sufficiently large finite domain and then diagonalized the 
obtained matrices. The solution of the ODEs and the diagonalization were per- 
formed using Mathematica 3.0. In all the cases the accuracy of the numerical 
solutions is better than a few per cent. 

6. 1 Properties of the static spike autosoliton 

In Sec. 3.1 we found that the static spike ASs can form in the Gray-Scott 
model when e <C 1. In one dimension these AS are stable when a > a = 0.33 
in the whole region of its existence A b < A < A d . Our numerical simulations 
of Eqs. (2.15) and (2.16) show that the value of A b is given by Eq. (3.6) with 
very good accuracy (less than 1%) already for e < 0.1, while the value of 
Ad approaches its asymptotic value Ad = 1.35 for e < 0.01 and increases 
somewhat for larger e (for example, for e = 0.05 we found A d = 1.48). These 
results are robust against decreasing e from 0.01 to 0.001 and smaller. For 
A b < A < 1 the distributions of 6 and r\ are described by Eqs. (3.3) and (3.7) 
with accuracy better than 10% for e < 0.1. 

The numerical solution of Eqs. (2.15) and (2.16) for e = 0.05 and A = 1.4 in 
the form of the one-dimensional static spike AS is presented in Fig. 23. The 
shape of the distributions of 9 and rj remains qualitatively the same all the 
way up to Ad = 1.48, but when A approaches A d the spike widens somewhat 
and its maximum becomes flatter. One can see that the potential V from 
the nonlinear eigenvalue problem [see Eq. (3.18)] shown in Fig. 23(b) indeed 
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assumes a complicated shape predicted by the asymptotic theory (Sec. 3). 
According to the asymptotic theory, the value of 6> max should decrease with the 
increase of A for A > 1.2 (see Fig. 5), and at A > A d the solution in the form 
of a single spike should disappear. This can also be seen from the simulations 
of the one-dimensional Gray-Scott model with e = 0.05 (Fig. 24). Observe the 
similarity between Fig. 24 and Fig. 5 obtained from the asymptotic theory. 
Also observe that the values of 6* max for A > 1 are smaller for finite values of e 
than those predicted by the asymptotic theory. When the value of A becomes 
greater than A d , a dip in the distribution of 9 appears in the center of the 
spike followed by a rapid decrease of the value of 9 there (local breakdown) 
resulting in the AS splitting into two spikes. These spikes go away from each 
other and then split again, until the system becomes filled with a periodic 
pattern of spikes. This process of self-replication of the static spike ASs in the 
Gray-Scott model was studied in detail in [50,55]. 

From the physical point of view, the existence of the static AS in which for 
e = l/L <C 1 the distribution of 9 = X/X has the form of a spike [see 
Fig. 23(a)] is determined by the diffusion processes and the special features 
of the autocatalytic reaction in Eq. (2.1) (see Sec. 2). Due to this reaction, 
the self-production of substance X occurs accompanied by the decrease of the 
amount of substance Y . The rate of this reaction is R ~ X 2 Y, so seemingly 
the concentration of substance Y should rapidly decrease. This, however, does 
not occur, and the high rate of the reaction in the spike is maintained by 
the strong diffusion influx of substance Y from the neighboring regions of size 
of order L. This ensures a finite concentration of substance Y in the spike 
which is necessary to maintain a high rate of the self-production of substance 
X. The diffusion current of substance Y is roughly Y Q /L, that is, it increases 
with A ~ Y [see Eq. (2.10)], but the reaction rate R grows as X 2 Y, so due 
to the increase of R with the increase of X ~ A ~ Y the concentration of 
Y ~ A~ 2 ~ Yq 2 [see Eqs. (3.5)] decreases faster in the spike as the value of 
Yq is increased. Therefore, for large enough Y the diffusion current cannot 
maintain high reaction rate in the spike. This explains why for large enough 
A ~ Y the value of 9 ~ X starts decreasing with the increase of A and at 
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Fig. 25. Dependences 6 , max (t) and w m m(t) for the static spike AS at the onset of the 
pulsations. Results of the numerical solution of Eqs. (2.15) and (2.16) for e = 0.05, 
a = 0.024 and A = 1.2 

some A = A d the solution in the form of the static spike AS disappears. 

6.2 Formation and collapse of the pulsating autosoliton 

In Sec. 5 we found that for a < 0.33 the static spike AS in the one-dimensional 
Gray-Scott model may become unstable with respect to the pulsations, with 
the critical value of a rapidly decreasing with the increase of A [see Eq. (5.12)]. 
Our numerical simulations showed that in a sufficiently large system this in- 
stability leads to the growth of self-oscillations of the AS amplitude. In other 
words, this instability leads to the transformation of the static AS into pulsat- 
ing [9-11]. Figure 25 presents the results of the simulation of such a process 
for e = 0.05, a = 0.024, and A = 1.2. The value of ct w = 0.0255 found in 
the simulations for these parameters agrees with the results of Sec. 5 with 
accuracy better than 5% (see Fig. 19). Also, the critical frequency agrees with 
the results of Sec. 5 within 1%. Note that for these values of A the value of 
q u given by Eq. (5.12) (which is also predicted in [51]) turns out to be off by 
a factor of 4. According to the results of the simulations, the pulsating AS 
seems to be unstable for all values of the parameters. After a few oscillations 
of its amplitude, it collapses into the homogeneous state (Fig. 25). 

6. 3 Properties of the traveling spike autosolitons 

The condition a C 1 means that the inhibitor Y is much slower than the 
activator X. The appearance of pulsation instability is due to the fact that 
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Fig. 26. Two types of traveling spike AS: case e » a 1 / 2 (a) and e ^ 
of the numerical solution of Eqs. (2.8) and (2.9). In (a) L = 0, a 
length is measured in the units of I. In (b) e = 0.05, a = 0.05, A 
measured in the units of L. 
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because of its sluggishness the inhibitor cannot control rapidly varying fluctu- 
ations of the activator. The lag in the inhibitor reaction also determines the 
existence of the traveling spike AS in the Gray-Scott model at a < 1 and 
e > a 1 / 2 . The supply of the substance Y into the spike of the traveling AS 
which is necessary to maintain high reaction rate R occurs by the AS running 
on the regions with high values of Y . The situation here resembles combustion 
process. The numerical simulations of Eqs. (2.15) and (2.16) in one dimension 
confirm the conclusions of Sec. 4 about the existence of the traveling spike 
ASs. Figure 26(a) shows the distributions of 6 and r\ in the form of an ultra- 
fast traveling spike AS for e = oo, a = 0.05, and A = 2. The speed of this 
AS was found to be c = 7.3, which agrees within 5% with that of Eq. (4.7). 
Also, the shape of the traveling spike AS is in agreement with the predictions 
of Sec. 4.1. 

When the value of e becomes smaller, or, more precisely, the value of L becomes 
sufficiently large, the concentration of substance Y in the front of the traveling 
spike decreases due to diffusion, so the speed of the AS decreases. This explains 
the possibility of the existence of two types of the traveling spike ASs studied 
in Sec. 4: the ultrafast traveling spike AS, which forms at e > a 1 / 2 , and 
the slower traveling spike AS, which forms at e 2 < a < e. The numerical 
simulations confirm this statement. Figure 26(b) shows the distributions of 6 
and rj in the form of a slower traveling spike AS obtained from the numerical 
solution of Eqs. (2.15) and (2.16) for e = 0.05, a = 0.05 and A = 2. The 
speed of this AS was found to be c = 1.26, much smaller than the speed of the 
ultrafast traveling AS discussed in the preceding paragraph. The shape of the 
slower traveling spike AS agrees with that found in the asymptotic theory (see 
Sec. 4.2). Figure 27 shows the dependence of the AS speed c on A at different 
values of a obtained from the numerical solution of Eqs. (2.15) and (2.16) 
with e = 0.05. Note that the curves c(A) in Fig. 27 terminate at sufficiently 
large values of A. This is due to the fact that when the value of A exceeds 
a certain critical value which depends on a, the traveling AS starts splitting 
as it moves. This is in agreement with the predictions of Sec. 4.2.3 about the 
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Fig. 27. Dependences c(A) for the traveling spike AS at different values of a obtained 
from the numerical solution of Eqs. (2.15) and (2.16) with e = 0.05. 

disappearance of the solution in the form of the traveling AS due to the onset 
of the oscillations in the back of the spike for sufficiently large values of A. 

Two ultrafast traveling ASs moving towards each other annihilate, what fol- 
lows from the physics of their existence. A much more diverse situation is 
realized in the case e 2 < a < e when the slower traveling spike AS exist. 
Here the ASs moving towards each other can annihilate before colliding or 
bounce off each other and start traveling in the opposite direction as a re- 
sult of the interaction via the diffusion of the inhibitor (a diffusion precursor). 
Also, as was shown in Sec. 5.1.7, when A > 1.58e 1//6 , the static spike AS may 
spontaneously transform into traveling when the value of a is decreased. This 
phenomenon was studied by Osipov and Severtsev in [45] in a simplified ver- 
sion of the Gray-Scott model and is observed in our simulations of Eqs. (2.15) 
and (2.16) as well. Comparing Fig. 27 with Fig. 24, one can see that there 
exists a parameter region where it is possible to excite both the static and the 
traveling ASs simultaneously. 

According to Fig. 27, when the value of A becomes sufficiently close to A d) the 
speed of the traveling spike AS may go to zero for a range of a. When e = 0.05, 
this happens when 0.03 < a < 0.06. According to the simulations, at e = 0.05 
and A = 1.34 the bifurcation of the static and the traveling ASs changes from 
subcritical to supercritical, so at higher values of A their coexistence is no 
longer possible. 

In our numerical simulations we also found the phenomenon of self-replication 
of the traveling spike ASs at large enough values of A. This phenomenon is 
illustrated in Fig. 28 which shows the density plots of the distributions of 9 
and i] as functions of x (horizontal axis) and t (vertical axis) for e = 0.05, 
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A = 2, and several values of a. One can see that for a = 0.05 the traveling 
AS goes back and forth between the system boundaries bouncing elastically 
off of them. When the value of a is increased to a = 0.07, a collision with 
the boundary stimulates a splitting event, which is followed by consecutive 
splitting of the newborn traveling spikes. Eventually, the system becomes filled 
with a stationary periodic pattern of spikes. Note that this pattern is stable in 
spite of the small value of a and the large value of A at which the solution in 
the form of the static spike AS does not exist. When the value of a is increased 
to a — 0.075, the traveling AS starts splitting as it moves, and when a = 0.1, 
splitting occurs more frequently. Observe that a small increase in the value 
of a leads to a significant decrease of the distance between the consecutive 
splitting events. Also, observe that despite splitting, the speed of the leading 
traveling spike remains practically constant. This is in agreement with the 
results of Sec. 4.2.3 which attribute the splitting to the onset of the oscillatory 
behavior in the tail of the traveling spike AS. The behavior of 9 and rj in 
the back of the spike only weakly affects the front of the spike, whose shape 
determines the speed of the traveling spike AS for sufficiently small values of 
a and sufficiently large values of A (see Sec. 4.2.2). Notice that self-replication 
of the traveling spike AS and wave reflection in the Gray-Scott model was first 
observed numerically in [56]. 

The results on the existence of the static and the traveling ASs, and the 
stability of the static AS in one dimension obtained in the previous sections, 
together with the results of the numerical simulations, are summarized in Fig. 
29. This figure shows the domains of existence and the instability lines for the 
ASs in the In a — In A plane in the limit e — > 0. 



7 Pattern formations scenarios in two dimensions 

Up to now, we studied the dynamics of the spike patterns in the one-dimensional 
Gray-Scott model. Now we are going to study the pattern formation scenarios 
in the two-dimensional Gray-Scott model. 

7.1 Granulation of the one- dimensional static spike autosoliton 

In Sec. 5.1.4 we showed that in two dimensions the static spike AS in the 
form of a stripe is unstable in the whole region of its existence A b < A < A d 
with respect to the corrugation instability, that is, the fluctuation 59q with 
k ~ 0.6Z -1 , even for a > 0.33. The growth of such a short-wave fluctuation 
should lead to the granulation of the stripe into small spots of size of order I. It 
is obvious that any wriggled stripe will also granulate into small spots. Figure 
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Fig. 28. Self-replication of the traveling spike AS. Results of the numerical solution 
of Eqs. (2.15) and (2.16) with e = 0.05, A = 2; a = 0.05 (a), a = 0.07 (b), a = 0.075 
(c), a = 0.1 (d). Darker areas correspond to lower values of 6 and rj. 
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Fig. 29. Diagram of the existence and the stability of different ASs in one dimension 
in the limit e — » 0. 
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Fig. 30. Granulation of a stripe. Results of the numerical solution of Eqs. (2.15) and 
(2.16) with e = 0.05, a = 0.5, and A = 2. The system is 2.5 x 2.5. The shades of 
gray show the distribution of rj. The spots show the regions where 9 > 10. 



30 shows such a process obtained from the simulations with e = 0.05, a = 0.5 
and A = 2. One can see that the stripe indeed granulates to spots of size 
of order I which then go away from each other until they become uniformly 
distributed across the system. Self-replication of spots may occur during this 
process (see below). 



7.2 Properities of the radially-symmetric static spike autosoliton 



As was shown in Sec. 3.3 and Sec. 5.3, static radially-symmetric AS in two 
dimensions exists in a relatively narrow range of the values of A ~ e. From our 
numerical simulations we see that at e = 0.05 and sufficiently large values of a 
a localized stimulus applied to the system at t = evolves into a stable static 
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Fig. 31. Self-replication of spots in two dimensions. Results of the numerical solution 
of Eqs. (2.15) and (2.16) with e = 0.05, a = 0.5, and A = 2. The system is 5 x 5. 
The shades of gray show the distribution of m The spots show the regions where 
9 > 10. 

radially-symmetric spike AS at 0.38 < A < 0.65. When a becomes sufficiently 
small, this region becomes even narrower since the AS becomes unstable with 
respect to the pulsations at sufficiently small A. If this is the case, an initial 
stimulus may produce a state which is close to a radially-symmetric AS, which 
after several pulsations will collapse. The process here is similar to the onset of 
the pulsations of the one- dimensional AS shown in Fig. 25. When the value of 
a becomes smaller than some value a c , static radially-symmetric AS becomes 
unstable for all values of A and can no longer be excited. 

7.3 Self-replication of the static radially-symmetric autosoliton 

As was shown in Sec. 5.2 and Sec. 5.3, when the value of A becomes greater 
than some critical value A&, the radially-symmetric static spike AS looses 
stability with respect to the radially non-symmetric fluctuations. The growth 
of such fluctuations leads to splitting and self-replication of the AS. In our 
simulations we found that for e = 0.05 and sufficiently large a static radially- 
symmetric spike AS becomes unstable and self-replicates at A > 0.7. Such a 
process for e = 0.05, a = 0.5, and A = 2 is shown in Fig. 31. From this figure 
one can see that the initial condition in the form of a rectangle of size of a few 
I splits into four (which is due to the rectangular shape of the initial condition 
and the fact that the value of A is well above A c2 ), and then the newborn 
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Fig. 32. Distribution of 9 in the simulation of Fig. 31 at t = 30. 



spots go on splitting until the system gets filled with an irregular arrangement 
of spots, with the characteristic distance between the spots much less than L. 
We would like to emphasize that the patterns observed in our simulations are 
essentially different from the domain patterns that form in N-systems [20-24] . 
The distributions of the activator in our simulations consist of the small spots 
instead of the sharp interfaces, and in the spots they are close to those in the 
radially-symmetric static spike AS. This is illustrated in Fig. 32 which shows 
the distribution of 6 at one of the moment of the simulation shown in Fig. 31. 

When the value of a becomes of order e, the dynamics of splitting significantly 
changes. Figure 33 shows the evolution of the system with e = 0.05, a = 0.1, 
A = 3, and a localized initial condition. As was already mentioned in Sec. 
6, a decrease in a means greater sluggishness of the inhibitor, so the pieces 
that form after splitting of an initial spot can go a greater distance apart and 
become more elongated than in Fig. 31 (where a ^> e). The state that forms 
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Fig. 33. Splitting as a result of the formation and the breakdown of a quasi 
one-dimensional wave. Results of the numerical solution of Eqs. (2.15) and (2.16) 
with e = 0.05, a = 0.1, and A = 3. The system is 5 x 5. The shades of gray show 
the distribution of r\. The spots show the regions where 9 > 10. 

here is close to a torn- up quasi one-dimensional wave of width of order I. This 
is natural to expect since, as we showed in Sec. 4.2 and Sec. 6, the traveling 
spike ASs are realized in the system for these values of the parameters. In this 
case the formation of new spots occurs as a result of their pinching off of the 
tips of the quasi one- dimensional wave pieces, that is, they sort of drip off from 
them. This process is illustrated in Fig. 34 which shows the distribution of 9 
at one moment of the simulation in Fig. 33. The spots that drip off from the 
wave pieces can further transform into quasi one-dimensional waves which in 
turn break up. As a result, the system becomes filled with a stable stationary 
pattern of spots, just as in the case of large a (see the last in Fig. 33). 

7.4 Spatio-temporal chaos 

For small enough values of a and A we were able to observe spatio-temporal 
chaos. Figure 35 shows the development of a chaotic pattern at e = 0.1, 
a = 0.04, and A — 1. This pattern does not transform to a stationary pattern 
of spots even for very long simulation times (t > 100). The stochastization 
of the pattern is caused by random splitting of spots and the disappearance 
of some of the spots due to their annihilation upon collision with the bigger 
spots. We would expect that these effects will be most pronounced at a ~ e 2 
and A ~ A c2 when the static radially-symmetric AS is close to the instabilities 
with respect to the pulsations and the onset of the traveling motion (Fig. 22) 
and is unstable with respect to splitting. Notice that such chaotic patterns are 
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Fig. 34. Distribution of 9 in the simulation of Fig. 33 at t = 2.44. 

also observed in semiconductor structures [12], combustion systems [14], and 
chemical systems [15]. 



7.5 Radially diverging waves 

For larger values of A and sufficiently small a an initially localized spot trans- 
forms into a circular stripe of growing radius (Fig. 36). When a is very small 
(a < e 2 ), this wave does not tear up and disappears at the system boundary 
[Fig. 36(a)]. When, on the other hand, we have a ~ e, the wave rebounds 
from the boundary and parts of it annihilate, so at some moment it tears 
up [Fig. 36(b)]. The tips of this wave get repelled from the boundary, so the 
wave propagates across the system until it annihilates upon collision with the 
boundary at its opposite side. 
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Fig. 35. Formation of spatio-temporal chaos. Results of the numerical solution of 
Eqs. (2.15) and (2.16) with e = 0.1, a = 0.04, and A = 1. The system is 10 x 10. 
The shades of gray show the distribution of rj. The spots show the regions where 
9 > 5. 



7.6 Spike spiral wave 



As was shown in Sec. 4.1 and Sec. 6, the ultrafast traveling spike AS can 
form in one dimension when a < 1 and e = oo (or L — 0). In our numerical 
simulations we found that in two dimensions the solution in the form of the 
traveling stripe is stable in a wide range of A, so it is natural to expect 
that at the same parameters it is possible to excite a steadily rotating spiral 
wave. The formation of such a wave at a = 0.1 and A = 2 is shown in Fig. 
37. We would like to emphasize that this wave is essentially different from 
the spiral waves observed in N-systems in that in the cross-section it has 
a form of a narrow spike and does not have a front and a back separated 
by a large distance, as is the case in N-systems (see, for example, [2,3,5-7]). 
Note that the spiral wave shown in Fig. 37 is precisely the kind of the wave 
that is observed in the experiments on Belousov-Zhabotinsky reaction, and 
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Fig. 36. Propagation of waves in the Gray-Scott model. Results of the numerical 
solution of Eqs. (2.15) and (2.16) with (a) e = 0.2, a = 0.05, and A = 1.5; (b) 
e = 0.1, a = 0.05, A = 1.5. In (a) the system is 10 x 10; in (b) the system is 5 x 5. 
The shades of gray show the distribution of ij. The spots show the regions where 
9 > 10. 

in the cardiac tissue, where such waves are thought to cause abnormal heart 
rhythms (see, for example, [3,5] and references therein). Recently, we proposed 
an asymptotic theory of these spiral wave in the Gray-Scott model [57]. 



8 Discussion 



While the work on this paper was being done, a number of publications on the 
Gray-Scott model appeared in the literature. The publications that are most 
relevant to our analysis are those of [50,51]. 

In [51] Doelman et al. present an asymptotic study of the static AS and peri- 
odic strata in the one- dimensional Gray-Scott model. They perform a singular 
perturbation analysis of the localized and spatially periodic stationary solu- 
tions in a limited region of the parameter space. Their results are in agreement 
with ours in this parameter region. However, the conclusions of the authors 
may be somewhat confusing since they do not use the natural scaling of the 
Gray-Scott model given by Eqs. (2.15) and (2.16). Instead, they introduce 
such dimensionless quantities that the original Eqs. (2.3) and (2.4) become 

4 = l6l-UV* + 5>a(l-U), (8.1) 
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Fig. 37. Formation of a steadily rotating spiral wave. Distribution of obtained 
from the numerical solution of Eqs. (2.8) and (2.9) with L = 0, a = 0.1, and A = 2. 
Length and time are measured in the units of / and re, respectively. The system is 
100 x 100. 



£ = <£ + <7V-^V, (8.2) 

where U is the inhibitor, V is the activator, S is the small parameter, a and b 
are constants of order one, and (3 (— ^ from the first of [51]) is a parameter 
that lies in the interval from zero to one. The solutions are studied in the limit 
5^0. 

It is clear that the scaling introduced in Eqs. (8.1) and (8.2) does not represent 
the true time and length scales of the problem, which are represented by Eqs. 
(2.15) and (2.16) of our paper. It is easy to see that Eqs. (8.1) and (8.2) can 
be reduced to Eqs. (2.15) and (2.16) with 



e = ^\/?, « = 5 2 "^, A = 6^. (8.3) 
V b b b 
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Observe that the scaling in Eqs. (8.1) and (8.2) forces certain relations be- 
tween the true parameters e, a, and A, so in fact such a choice of scaling 
significantly restricts the parameter space studied by Doelman et al, leading 
them to incorrect conclusions about non-existence of certain types of solutions. 

According to the results of the first paper in [51], the static one-dimensional 
AS exists only when (3 G [0, 1) as 5 — > 0. This statement is based on the fact 
that the coefficient in front of the last term in the right-hand side of Eq. (8.1) 
scales as 5 2 . On the other hand, as we showed in Sec. 3.1, this AS in fact exists 
in a wider range of the parameters as long as e -C 1, so in fact such a scaling 
is not a necessary condition for the AS existence. This fact was noticed by the 
authors in the more recent paper in [51]. 

The case (3 = is equivalent to A ~ Ab in our notation. Our result that the 
solution exists only at A > Ab = \/12e is in agreement with the result of [51] 
that the solution exists at a > 1446 3 . Notice that the analogous method of 
finding Ab for the ASs in systems of small size was introduced two decades 
ago by Kerner and Osipov [16]. There they performed an analysis of a very 
similar model — the Brusselator. The results obtained by them at A ~ A& 
for the Brusselator differ from those presented in our paper only by numeri- 
cal coefficients. A calculation similar to that of Sec. 3.1.1 was performed by 
Dubitskii, Kerner, and Osipov for the Gierer-Meinhardt model [42] (see also 
[11]). Osipov and Severtsev also performed such an analysis for a simplified 
version of the Gray-Scott model, which are also very similar to the results of 
the first part of Sec. 3.1 [45]. The advantage of the method of Doelman et al., 
however, is that it rigorously shows the existence of the static spike AS for 
A b < A < 1 and e < 1. 

Doelman et al. claim that the solution obtained at (3 = is valid with accuracy 
to 8 ~ e 1 / 2 . On the other hand, the analysis of Sec. 3.1 suggests that this 
solution is actually valid with a much better accuracy e. In fact, e is the true 
small parameter of the problem, and the asymptotic procedure for obtaining 
the solutions for e 1 can be constructed regardless of all other parameters 
of the problem with accuracy better than e 1 / 2 (see the discussion at the end 
of Sec. 3). Note that this procedure for general V- and A-systems was first 
introduced by Dubitskii, Kerner, and Osipov in [42] (see also [11]). 

For (3 > Doelman et al. claim that the solution obtained by them is accurate 
to order <5 1_/3 . They forget that although the accuracy with which the asymp- 
totic equations describe the actual distribution of V increases as (3 decreases, 
the accuracy of the matching condition, namely, the fact that the true U(0) 
is non-zero, decreases as (3 decreases.^ It is not difficult to see that the latter 



5 In the first of [51] this is equivalent to the approximation of Eq. (3.6) by Eq. (3.7) 
of that paper. See also the remark after Theorem 4.3 there. 
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actually gives corrections of order 5 2 which become of order 1 as (3 — > 0, 
as should be expected. This means that the best accuracy with which the 

3 

asymptotic procedure of Doelman et al. is valid is 5~>, what in our notation 
corresponds to e 1 / 3 . On the other hand, as can be seen from Sec. 3.1, the 
asymptotic procedure introduced by us always gives an accuracy better than 
e 1//2 , and in the parameter regions of interest the accuracy is of order e. Fur- 
thermore, our asymptotic procedure remains valid even for A ~ 1, which is 
equivalent to f3 — 1 of [51], where the approach of the latter fails completely, 
thus being able to quantitatively describe the local breakdown and the dis- 
appearance of the solution at A = A d . Notice that in the second of [51] the 
authors used a different approach to qualitatively analyze the disappearance 
of solution at A > 1. 

Doelman et al. make a statement that the traveling ASs do not exist in the 
Gray-Scott model. Although this statement is true for the scaling in Eqs. 
(8.1) and (8.2), it is generally wrong. In fact, as we showed in Sec. 4, there 
are two qualitatively different traveling spike ASs that are realized in the 
Gray-Scott model. The ultrafast traveling spike AS are realized when a < e 2 
and a 1 / 2 < A < a -1 / 2 , and the slower traveling spike AS are realized when 
g2 ~ a ~ e an d ae_1 ~ A < cT x l 2 . These parameter regimes lie outside of 
those studied in [51]. The more so, we find coexistence of the traveling spike 
AS with static. 

In the second of [51] Doelman et al. study the stability of the one-dimensional 
static spike AS. Their analysis is based on the equations of the type of Eq. 
(5.9) and (5.11). Equations of this type were studied by Kerner and Osipov 
in the context of the stability of the spike AS in the systems of small size [9- 
11,16]. As we showed in the appendices A and B, a straightforward extension 
of the method of Kerner and Osipov allows rigorous analytic studies of these 
equations. Doelman et al. found in the scaling of Eqs. (8.1) and (8.2) that 
the AS becomes unstable with respect to the pulsations when j3 — |, what 
corresponds to A u ~ e 2 / 7 . Their results agree with those obtained by us in 
Sec. 5.1.3 in the case A b <C A <C A d [see Eq. (5.12)]. However, as can be seen 
from Sec. 5.1.3 and 6 of our paper, in practice this result has a very limited 
applicability, so one should use the results presented in Fig. 19 of our paper 
to find the values of and u . 

Moreover, this is not the only possibility for the instability of the static spike 
AS for the general choice of the parameters. What we showed in Sec. 5.1 is that, 
first, the instability with respect to the pulsations is realized in a wider range 
of the system parameters, and, second, there is another instability which leads 
to the transformation of the static spike AS into traveling, which for A > e 1 / 6 
precedes the pulsation instability. In general, we performed a complete analysis 
of the stability of the static spike AS in one dimension. 
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Reynolds, Ponce-Dawson, and Pearson derived asymptotic equation of motion 
for the spike AS in the Gray-Scott model with A ~ 1, a ~ 1, and e 1 (in 
our notation). The reason the AS moves is its interaction with the boundary. 
However, for A > A d the peculiarities of the internal dynamics result in the 
breakdown of the asymptotic description which the authors attribute to split- 
ting and self-replication of the AS. The equations obtained by the authors for 
the inner region in fact coincide with Eqs. (4.19) and (4.20) obtained by us in 
Sec. 4.2 for the traveling spike AS in the case e 2 < a < e. They also found the 
value of A d which agrees with the one obtained by us. Thus, the peculiarities 
noticed in [50] should be important in the case of the traveling spike AS in 
the system of infinite extent and should explain the dynamics of splitting of 
the traveling AS for A ~ 1 and a ~ e. 

Notice that both the authors of [50] and the authors of [51] studied self- 
replication and splitting of the static spike ASs in one dimension and their re- 
sults cannot be used to explain self-replication in real chemical systems which 
are higher- dimensional and in which the radially-symmetric static spike ASs 
(spots) rather than the one-dimensional static spike ASs (pulses) form. Indeed, 
as we showed in the Sec. 5.1, 5.2, 6 and 7, in one dimension splitting occurs as 
a result of the local breakdown in the AS center, whereas in higher dimensions 
splitting is the result of the buildup of radially non-symmetric fluctuations and 
occurs at different values of the parameters. Moreover, the local breakdown 
in the case of the higher-dimensional radially-symmetric static spike AS can- 
not be realized at all. Note that the same situation is realized in a different 
class of reaction-diffusion systems studied by us (N-systems) where splitting 
is also shown to be the consequence of the buildup of radially non-symmetric 
fluctuations [21,22,40]. Thus, for e 1 self-replication and splitting in both 
N-systems and the Gray-Scott model are driven by qualitatively the same 
mechanism. 

Let us comment on the relationship between the numerical simulations of the 
two-dimensional Gray-Scott model performed by Pearson in [26] and those 
of Sec. 7 performed by us. Pearson uses a different non-dimensionalization 
of the Gray-Scott model, which has the following correspondence with our 
parameters: e 2 = ^rf^y, a = ^rj^, A = (see [26]). It is not difficult to 
see that for the simulations of Pearson e ~ 0.45, a ~ 0.4, and A ~ 2, so his 
choice of the parameters corresponds to e ~ 1 and a ~ 1. This is different 
from our simulations for which e < 1 and/or a C 1. Note that the stationary 
patterns in the Gray-Scott model with e ~ 1 should actually resemble those 
forming in N-systems [22]. Indeed, if one introduces the new variables 9 = -| 

2 

and 77 = 77 + ^9, after simple algebra one can write Eqs. (2.15) and (2.16) as 



a^- = e 2 A9 + A 2 9 2 f) - e 2 A 2 9 3 - 9 (8.4) 
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From Eqs. (8.4) and (8.5) one can see that for e ~ 1 the nullcline of the 
equation for 6 is actually N-like, and the coupling between 6 and rj becomes 
linear, so the stationary patterns in this case should in fact look like those 
forming in N-systems [22]. This is the reason why Pearson observed the sta- 
tionary labyrinthine patterns, while in our simulations any stripe-like pattern 
always granulates into spots. Also, the collective oscillations of the space-filling 
patterns observed by Pearson are similar to the collective oscillations of the 
domain patterns in N-systems [41]. Pearson did not see the spiral waves be- 
cause in his simulations a ~ e 2 and the spirals break up as they form. We do 
not see any phase turbulence since in our simulations the system is far away 
from the Hopf bifurcation of the homogeneous state 6^3, ^3. The rest of the 
patterns observed in [26] are similar to those observed by us. 

Finally, we will mention recent studies by Hale, Peletier, and Troy [58]. They 
found the exact solutions in the form of the solitary pulses and fronts and 
analyzed their stability in the Gray-Scott model with e = 1 and a = e 2 . 
The reason this can be done exactly for these values of e and a is because 
the equation for the activator and the inhibitor effectively decouple from each 
other [see Eq. (8.5)], so the model behaves as a scalar reaction-diffusion system 
[5,10,11]. It is then not surprising that the behavior of the system remains the 
same if e is close to, but not exactly equal to 1 [58]. 



9 Conclusion 

Let us now summarize the results of our analysis of the patterns in the Gray- 
Scott model. As was emphasized in Sec. 2, a unique feature of the Gray-Scott 
model is the fact that in it the homogeneous state Oh — 0, rjh — 1 is stable for 
all values of the parameters. However, in such a stable homogeneous system it 
is possible to excite various steady inhomogeneous states, including the self- 
sustained solitary inhomogeneous states, autosolitons (ASs), by applying a 
sufficiently strong external stimulus. The formation of these inhomogeneous 
states is due to the self-production of substance X (which plays the role of 
the activator) controlled by the other substance Y (which plays the role of 
the inhibitor). The properties of the patterns are determined by only three 
parameters: e, a, and A. The parameters e = l/L and a = Tg/r^ are the 
ratios of the characteristic length and time scales of the activator and the 
inhibitor, respectively, and the control parameter A determines the degree of 
the deviation of the system from equilibrium since it is proportional to the 
rate of supply of substance Y, which plays the role of "fuel" for the reaction in 
Eq. (2.1). We emphasize that for the same values of the system's parameters 
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Table 1 



Scaling of the main parameters of different types of ASs in the Gray-Scott model. 
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it is possible to excite different patterns by choosing the form of the stimulus, 
which will be stable in certain ranges of the parameters e, a, and A. At the 
stability margin the patterns spontaneously disappear or transform into the 
patterns of different kind. 



As follows from the general qualitative theory of the patterns in reaction- 
diffusion systems, the necessary condition for the existence of the persistent 
patterns of any kind is the smallness of the parameters e and/or a [9-11]. In 
this paper we took advantage of this fact and performed asymptotic analysis 
of the simplest patterns (ASs) in the limit when either of these parameters 
goes to zero. This analysis gave us the dependence of the main parameters of 
the ASs on the system's parameters and their ranges of existence and stability. 
Their asymptotic behavior is summarized in Table 1. The properties of the 
forming patterns strongly depend on the parameters of the system. Depending 
on the values of e and a one can distinguish three different cases. 

The first case corresponds to e < 1 and a > 1. As was expected from the 
general qualitative theory [9-11], for these values of e and a one can excite only 
the static spike ASs (Fig. 23). In one dimension these ASs are stable in a wide 
range of A from A^ ~ e 1 / 2 to A^ ~ 1 (Sec. 3.1). The characteristic size of the 
spike of the AS is determined by the diffusion length / of the activator and is 
practically independent of A. When the value of A is increased, the amplitude 
of the static spike AS grows from # max ~ e™ 1 / 2 for A ~ A b to # max ~ e _1 for 
A ~ 1. According to the general qualitative theory [9-11], at A = A b and 
A = Ad we must have d9 max /dA = oo. For this reason at A = A b the static 
spike AS, having a large amplitude, abruptly disappears. When A approaches 
Ad, the AS widens somewhat, until at A = Ad a local breakdown occurs in its 
center leading to the consecutive splitting of the static spike AS. As a result 
of this self-replication effect a static pattern consisting of a periodic array of 
spikes forms [50]. 

When the value of a is decreased below a ~ 1, the range of existence of the 
static spike AS narrows. When A is decreased, the AS looses its stability before 
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reaching the point A = A\> (Sec. 5). The instability is realized with respect 
to the pulsations leading to the AS collapse (Fig. 25). On the other hand, 
when A is increased, the static spike AS may destabilize and spontaneously 
transform into traveling before reaching the point A = A d (see also [45]). 

In the two- and the three-dimensional Gray-Scott model with eC 1 and a > 1 
one can excite the radially-symmetric static spike ASs of size of order / (Sec. 
3.2 and Sec. 3.3). The range of the values of A for which these ASs exist is 
very narrow for e 1 (see Table 1). At the point A = A b ~ e the static 
radially-symmetric spike AS, having large amplitude # max ~ e" 1 , abruptly 
disappears. The range of A at which the radially-symmetric static spike AS 
exist becomes even narrower for a ~ e 2 when the AS becomes unstable with 
respect to the pulsations (Sec. 5.2 and Sec. 5.3). On the other hand, when the 
value of A is increased to A = A c2 ~ e, the static radially-symmetric spike AS 
looses stability with respect to the radially non-symmetric fluctuations. As a 
result of the development of such fluctuations the AS splits into two, which 
then split in turn (self-replicate) until the system gets filled with a multispot 
pattern (Fig. 31). We would like to emphasize that for e 1 a spot (a state 
close to the radially-symmetric AS) is the dominant morphology, so that any 
localized initial state such as a stripe or a square first granulates into spots, and 
the evolution of the system is then governed by self-replication of these spots 
(Figs. 30, 31). Let us note that in N-systems one sees complex patterns in the 
form of wriggling stripes, connected and disconnected labyrinthine patterns 
[20-24], which are also observed in chemical experiments [15]. These patterns 
do not form from a localized stimulus in the Gray-Scott model with e < 1. 
Note, however, that when e ~ 1, the Gray-Scott model starts behaving like 
an N-system (see the end of Sec. 8) , so for these values of e such patterns can 
in fact be excited [26]. 

In the second case we have a C 1 and e 3> a 1 / 2 . In this case one can excite 
different kinds of self-sustained waves (autowaves) which have the form of 
the narrow spikes of size roughly I and the amplitude # max ~ a -1 . In one 
dimension the ultrafast traveling spike ASs are realized [Fig. 26(a)] whose 
speed c ~ AoT x l 2 x 1/tq, i. e., it is AaT 1 ! 2 times greater than the typical speed 
of the traveling AS in N-systems (Sec. 4.1). This ultrafast traveling spike AS 
can be excited in a wide range of A from A b ~ a 1 / 2 <C 1 to A d ~ cT 1 ! 2 ^> 1. 
In two dimensions, besides the ultrafast traveling spike AS, one can excite 
radially diverging waves (Fig. 36) and the steadily rotating spiral wave (Fig. 
37), which in the cross-section looks like the ultrafast traveling spike AS. The 
spiral wave observed by us is a new type of spiral waves, which is essentially 
different from those forming in N-systems. 

In the third case 6 < a € 1 the behavior of the patterns in the Gray-Scott 
model is most diverse. In one dimension, besides the static spike AS one can 
excite the traveling spike AS [Fig. 26(b)] whose speed decreases with the de- 
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crease of e. These ASs may elastically rebound upon collision with the bound- 
ary or with each other (Fig. 28). When a > e, they start self-replicating. The 
greater the value of a, the smaller the distance between the consecutive split- 
ting events in such a self-replication process. We would like to emphasize that 
as a result of this effect a static periodic pattern of spikes forms (Fig. 28). In 
two dimensions the traveling spike AS undergoes a transversal breakup lead- 
ing to the formation of the stationary multispot patterns (Fig. 33). Also, when 
the values of a and A are small enough, one can excite the turbulent patterns 
(Fig. 35). The chaotic behavior of the latter is due to the random creation of 
the new spots as a result of the self-replication and the annihilation of some of 
the spots as they collide with each other. This kind of turbulence is observed 
in chemical experiments [15] and is not unlike the one realized in N-systems 
[22]. 

Above we considered the properties of the spike ASs, which are the simplest 
patterns and are, therefore, the building blocks of the more complex patterns 
forming in the Gray-Scott model. Their properties allowed us to understand 
the pattern formation scenarios in the Gray-Scott model subjected to a lo- 
calized stimulus. We expect that the asymptotic methods developed by us in 
this paper can be successfully applied to other systems of this kind. Also, we 
found that in many instances a localized stimulus results in the formation of 
the complex patterns consisting of many strongly interacting AS-like states. 
The effect of interaction of such states in the complex patterns is an interesting 
open problem. We hope that the singular perturbation techniques developed 
in this paper will be useful for studying these interactions and their effects on 
the dynamics of the complex space-filling patterns in the Gray-Scott model 
and other similar models. 
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A Analysis of Eq. (5.9) 

Equation (5.9) is of the kind studied by Kerner and Osipov in the case of the 
ASs in systems of small size [9-11,16]. Here we perform a rigorous analysis of 
Eq. (5.9) using their method. 

Let us introduce the orthonormal basis set 59 n of the eigenfunctions of the 
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Schrddinger operator in the left-hand side of Eq. (5.9) Q 



d 2 { x 

+ 1 — 3 cosh - ' 



dx 2 



S9 n = \J9 n . (A.l) 



For the simplicity of notation we omitted the superscript in this section. As 
was already discussed in Sec. 5.1.1, this operator has three discrete eigenvalues 
[Eq. (5.7)] and a continuous spectrum for A n > 1. 

Assuming for a moment that the problem is considered on a large but finite 
domain, we can write the operators of Eq. (5.9) in this basis as 



B mn = (A„ - i)5 mn + C(l - j)b l m b r n , (A.2) 
where 5 mn is the Kronecker delta, 

+oo +oo 

b l n = J cosh -4 ^— J 59 n (x)dx, b r n = J 59 n (x)dx, (A.3) 

— oo — oo 

and 



Observe that C is a monotonically increasing function of A > A^. 
In terms of B mn Eq. (5.9) becomes 

det B mn = 0. (A.5) 

Note that since by symmetry b l n and b r n are identically zero for odd functions 
S9 n , we immediately conclude that these functions are the solutions of Eq. 
(5.9) with 7 n = X n corresponding to these functions. 

It is not difficult to show that because of the special form of the second matrix 
in Eq. (A.2) we have [9-11,16] 



There should be no confusion between the eigenfunctions 59 n of this section, which 
)rrespond t< 
of Eq. (5.6). 



correspond to 59 n °^ of Sec. 5.1 and 59 n of that section, which are the eigenfunctions 
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det B„ 



1 + C(1- 7 )E 



,,, A n 



7 



n(*»-7), 



(A.6) 



where a n = b l n b r n and the summation is over the even states only. Using Eq. 
(A.l), one can bring the expression for a n to a symmetric form which is con- 
venient for the further calculations 



2A r 



a n = 



A r , 



+00 

J cosh -2 (^j 56 n (x)dx 



(A.7) 



The values of a and a 2 can be calculated explicitly with the use of Eq. (5.7) 



a 



757T 2 

256 



a 2 



97T 2 

256' 



(A.8) 



The calculation of corresponding to the functions 56k of the continuous 
spectrum (with the wave vector k and A& = 1 + k 2 ) is rather involved. The 
functions 56k can be written as linear combinations of the real and the imag- 
inary parts of 



u(y) = (1 - y 2 ) ik F (2ik - 3, 2ik + 4, 2ifc + 1, ^) , (A.9) 

where y = tanh(x/2) and F(a, /3, 7, x) is the hypergeometric function [54], to 
obtain the even functions 56k- The functions 56k should be normalized in such 
a way that 56k(x) — > cos(A;a; ± 5) as x — > ±00. Then, after calculating the 
respective integrals, we arrive at 

_ 87r 2 fc 2 (fc 2 + 1) 

° fc ~ (16A; 4 + 40fc 2 + 9) sinh 2 (7rA;) 

Naturally, in the infinite domain one should replace the summation over the 
continuous spectrum in Eq. (A.6) by integration: J2 n ~~ * Jo°° V' 

To study the unstable solutions of Eq. (A. 5), we need to analyze the zeros of 
the function 

%) = l + C(l + i,)(-^ + -^+f , Qfc /f , ) (A.ll) 
K ' v ^Ao + icu A 2 + icj J 7r(l + k 2 + iu) J v ; 

in the lower half-plane of the complex frequency 00 = i'y. This can be done 
with the aid of the argument principle [9-11,16] which states that the number 



(A.10) 
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of zeros iV of the complex function D(iv) in this region of complex frequency 
iv is equal to 



N = P + —A arg£)(u;), (A.12) 

where P is the number of poles there and Aarg.D(u;) is the change of the 
argument of the function D(u) as iv winds around this region of the complex 
frequency counterclockwise. 

From the spectrum of the operator in Eq. (A.l) one can see that only the pole 
at iv — iAo lies in the lower-half plane of the complex frequency, so we have 
P = 1 [9-11,16]. Let us see how the function D(iv) with iv real varies as uo 
goes from +00 to —00. Since D{ui) is symmetric with respect to the real axis, 
one only needs to analyze the case of positive iv. At iv = 00 we have 



£>(oo) = 1 + C (a, + a 2 + J^-j > 0, (A.13) 

where we used the explicit expressions for a 0i 2,fc and A ,2 to calculate the sign 
of -D(oo) and evaluated the integral in this equation to be J °° n~ 1 a k dk ~ 0.12. 
On the other hand, at iv = we have 

l c<0 (A14) 



for A > At,. The latter expression can be obtained by recalling that at A = Aj,, 
for which C = 3/8, we have D(0) = (see Sec. 5.1.1), and C monotonically 
increases with A. 

It is not difficult to show that the imaginary part of D{iv) is negative for all 
iv > 0: 

Im D((v) a (\ — 1) 02(^2 — 1) 7 k 2 akdk 
Ctv \ 2 + lv 2 + \ 2 2 + u 2 + l tt[(1 + A; 2 ) 2 + ^ 2 ] 

dn(\n — 1) 1 / _x f k 2 akdk\ „ . . . 

< T2~7 2 + VT— - «2 A 2 - 1 + / — < 0. A.15 

The last inequality is obtained by using the explicit expressions for ao,2,fc ; Ao,2 
and the evaluation of the last integral J °° n~ 1 k 2 a k dk ~ 0.02. Note that because 
of the smallness of the contributions from the continuous spectrum one can 
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Fig. A.l. The behavior of the function D(uS) for the static AS (a) and for the small 
amplitude state corresponding to the "-" sign in Eq. (3.5). 

get very good approximations for the solutions of Eq. (5.9) by restricting S6 n 
to the discrete spectrum only (see also [45]). 

From all this we conclude that the function D(u) has the form shown in Fig. 
A. 1(a), so we have AaigD(u) = —2n for the static spike AS. This means 
that = and Eq. (5.9) does not have solutions with Re 7 < 0. Note that 
the same line of arguments shows that the asymptotic stability problem for 
the stationary solution with the smaller amplitude which corresponds to the 
"-" sign in Eq. (3.5) always has a solution with Re 7 < since in that case 
Aarg_D(oj) = [see Fig. A. 1(b)], so this small-amplitude solution is always 
unstable. These conclusions are in agreement with the general qualitative the- 
ory of the ASs [9-11]. 



B Analysis of Eq. (5.11) 

In this section we use the method of the previous section to analyze the solu- 
tions of Eq. (5.11). This method was used by Kerner and Osipov for studying 
the instabilities of the static ASs for small values of a in the systems of small 
size [9-11]. 

After introducing the orthonormal basis set of the eigenf unctions of Eq. (A.l), 
we get the following expression for B mn for Eq. (5.11): 
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Brnn — (A n + \Lo)8 mn + C ' —~ 

yw + a 



1 + b l ¥ 

u m u ni 



(B.l) 



where now 



C 



A 2 a 1 ' 2 
8e : 



(B.2) 



the rest is the same as in Eq. (A. 2), and a — > +0 [cf. Eq. (5.5), a will determine 
the proper winding direction, see below]. As in the previous section, we may 
write 



det Brnn. — 
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(B.3) 



where a n are given by Eq. (A. 7). To analyze the solutions of Eq. (A. 5) with 
this B mn , we will study the zeros of the function 



D(u) = 1 + C 
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dhdk 
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where a fc are given by Eq. (A. 10), in the lower half-plane of the complex 
frequency u by using the argument principle [Eq. (A. 12), in which, as before, 
P = 1]. Of course, as in the previous section, D{uj) should be symmetric with 
respect to the real axis. 

For uj > the real and the imaginary parts of D(u>) can be written as 



Re 



and 



2u 



+ + 



a A c 



C 

+oo(l - uj) 



\l + u 2 
a 



+ 



a 2 A 5 



\ 2 2 + u 2 
a 2 



+ 



(1 + k 2 )a k dk 



o 

oo 



\l +U 2 Al + CJ 2 



/vrffl 



tt[(1 + A; 2 ) 2 + cu 2 ] 
dkdk 



[(1 + k 2 ) 2 + LU 2 } 



:b.5) 



Im 



/2lu 



+{u-l) 



a 



a 2 



Ag + cu 2 Ai + cu 2 7 tt[(1 + A; 2 ) 2 + cu 2 ] 



dkdk 



aoAo 

a|w 



+ 



a 2 A 2 



X 2 + uj 2 



+ 



J 7T[fl 



(1 + k 2 )a k dk 



[(1 + A; 2 ) 2 + cj 2 ] 



(B.6) 
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Fig. B.l. Qualitative form of the function D(uj) for large values of C (a) and for 
small values of C (b). 

Using the explicit expressions for a 0;2 ,k and A ,2 from the previous section, it 
is not difficult to show that the expressions in the brackets above are negative 
for all values of uj. The analysis of Eq. (B.6) then shows that Im D(u) should 
change sign once when < uj < oo. Let us denote the value of uj at which this 
happens as uj . Note that according to Eq. (B.6), we must have cu < 1. 

From the definition of D(uj) one can see that 



oo 



D(u) -> 1 + C^=L \a + a 2 + f ^— ) , uj -»• ±oo. (B.7) 
■2\u\ V I * ) 

Since the expression in the bracket in this equation is positive, we will have 
Im D{uj) < for sufficiently large uj > 0. On the other hand, 



^^ Md^) , w _>± . (B.8) 
3V 2 I W I 

Therefore, for sufficiently small uj > we must have Im D(uo) > 0. Observe 
that Eq. (B.8) was obtained for a — 0. When a is small but finite, the two 
branches in Eq. (B.8) will actually get connected at Re D(uo) ~ — a^ 1 ^ 2 . Thus, 
the qualitative behavior of D(uo) should be the one shown in Fig. B.l Note 
that the function D(u) can be calculated numerically from Eq. (B.4) for any 
value of C and has indeed the form shown in Fig. (B.l). 

The number of zeros of D(u) in the lower half-plane of the complex frequency is 
determined by Re D(ojq). According to Eq. (B.4), if C is sufficiently small, the 
first term in Eq. (B.5) will dominate for uj = uj , so we will have Re D(uj ) > 0. 
In this case the change of the argument of D(uj) will be A axgD(uj) = 2ir [Fig. 
B.l (a)], so we will have N = 2 and therefore an instability. On the other hand, 
if C is large, we can neglect \f2ujjC at uj = uj in Eq. (B.5). From Eq. (B.5) 
and the fact that uj < 1 one can then see that Re D(uj ) < 0. In this case 
the change of the argument will be A aigD(uj) = —2n [Fig. B.l(b)], so the 
number of zeros in the lower half-plane is iV = 0, implying stability. From all 
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this we see that as the value of C is decreased, at some C = Co a complex- 
conjugate pair of unstable solutions of Eq. (5.11) appears signifying a Hopf 
bifurcation. The numerical analysis of Eq. (B.4) shows that C — 0.2837, what 
corresponds to a u ~ 5.15e 2 A~ 4 , and c^o — 0.534, in excellent agreement with 
the results of Sec. 5.1.3. 
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